计算几何模板（大神整理）

不知道出处，转载的转载。。。

计算几何模板

目录：

1.计算几何 2

1.1 注意 2

1.2几何公式 2

1.3 多边形 4

1.4多边形切割 7

1.5 浮点函数 8

1.6 面积 14

1.7球面 15

1.8三角形 18

1.9三维几何 21

1.10 凸包 29

水平序 29

极角序 30

卷包裹法 31

1.11 网格 33

1.12 圆 34

1.13 矢量运算求几何模板  36

1.14结构体表示几何图形 48

1.15四城部分几何模板 53

1.16  一些代码 55

1.16.1 最小圆覆盖_zju1450  55

1.16.2 直线旋转_两凸包的最短距离(poj3608) 59

1.16.3 扇形的重心 63

1.16.4 根据经度纬度求球面距离  64

1.16.5 多边形的重心 65

1.16.6 存不存在一个平面把两堆点分开(poj3643) 67

1.16.7 pku_3335_判断多边形的核是否存在 68

1.16.8 pku_2600_二分+圆的参数方程 75

1.16.9 pku_1151_矩形相交的面积 77

1.16.10 pku_1118_共线最多的点的个数 79

1.16.11 pku2826_线段围成的区域可储水量 81

1.16.12 Pick公式  85

1.16.13 N点中三个点组成三角形面积最大 87

1.16.14 直线关于圆的反射  90

1.16.15 pku2002_3432_N个点最多组成多少个正方形(hao) 95

1.16.16 pku1981_单位圆覆盖最多点(poj1981)CircleandPoints 98

1.16.17 pku3668_GameofLine_N个点最多确定多少互不平行的直线(poj3668) 100

1.16.18 求凸多边形直径  101

1.16.19 矩形面积并，周长并  103

1.16.20 pku2069 最小球覆盖  103

1.16.21 最大空凸包、最大空矩形  106

1.16.22 求圆和多边形的交  106

半平面交 110

Nlgn 110

N^2 112

1.计算几何

1.1 注意

1. 注意舍入方式(0.5的舍入方向);防止输出-0.

2. 几何题注意多测试不对称数据.

3. 整数几何注意xmult和dmult是否会出界;

   符点几何注意eps的使用.

4. 避免使用斜率;注意除数是否会为0.

5. 公式一定要化简后再代入.

6. 判断同一个2*PI域内两角度差应该是

   abs(a1-a2)<beta||abs(a1-a2)>pi+pi-beta;

   相等应该是

   abs(a1-a2)<eps||abs(a1-a2)>pi+pi-eps;

7. 需要的话尽量使用atan2,注意:atan2(0,0)=0,

   atan2(1,0)=pi/2,atan2(-1,0)=-pi/2,atan2(0,1)=0,atan2(0,-1)=pi.

8. cross product = |u|*|v|*sin(a)

   dot product = |u|*|v|*cos(a)

9. (P1-P0)x(P2-P0)结果的意义:

   正: <P0,P1>在<P0,P2>顺时针(0,pi)内

   负: <P0,P1>在<P0,P2>逆时针(0,pi)内

   0 : <P0,P1>,<P0,P2>共线,夹角为0或pi

10. 误差限缺省使用1e-8!

1.2几何公式

三角形:

1. 半周长 P=(a+b+c)/2

2. 面积 S=aHa/2=absin(C)/2=sqrt(P(P-a)(P-b)(P-c))

3. 中线 Ma=sqrt(2(b^2+c^2)-a^2)/2=sqrt(b^2+c^2+2bccos(A))/2

4. 角平分线 Ta=sqrt(bc((b+c)^2-a^2))/(b+c)=2bccos(A/2)/(b+c)

5. 高线 Ha=bsin(C)=csin(B)=sqrt(b^2-((a^2+b^2-c^2)/(2a))^2)

6. 内切圆半径 r=S/P=asin(B/2)sin(C/2)/sin((B+C)/2)

               =4Rsin(A/2)sin(B/2)sin(C/2)=sqrt((P-a)(P-b)(P-c)/P)

               =Ptan(A/2)tan(B/2)tan(C/2)

7. 外接圆半径 R=abc/(4S)=a/(2sin(A))=b/(2sin(B))=c/(2sin(C))

四边形:

D1,D2为对角线,M对角线中点连线,A为对角线夹角

1. a^2+b^2+c^2+d^2=D1^2+D2^2+4M^2

2. S=D1D2sin(A)/2

(以下对圆的内接四边形)

3. ac+bd=D1D2

4. S=sqrt((P-a)(P-b)(P-c)(P-d)),P为半周长

正n边形:

R为外接圆半径,r为内切圆半径

1. 中心角 A=2PI/n

2. 内角 C=(n-2)PI/n

3. 边长 a=2sqrt(R^2-r^2)=2Rsin(A/2)=2rtan(A/2)

4. 面积 S=nar/2=nr^2tan(A/2)=nR^2sin(A)/2=na^2/(4tan(A/2))

圆:

1. 弧长 l=rA

2. 弦长 a=2sqrt(2hr-h^2)=2rsin(A/2)

3. 弓形高 h=r-sqrt(r^2-a^2/4)=r(1-cos(A/2))=atan(A/4)/2

4. 扇形面积 S1=rl/2=r^2A/2

5. 弓形面积 S2=(rl-a(r-h))/2=r^2(A-sin(A))/2

棱柱:

1. 体积 V=Ah,A为底面积,h为高

2. 侧面积 S=lp,l为棱长,p为直截面周长

3. 全面积 T=S+2A

棱锥:

1. 体积 V=Ah/3,A为底面积,h为高

(以下对正棱锥)

2. 侧面积 S=lp/2,l为斜高,p为底面周长

3. 全面积 T=S+A

棱台:

1. 体积 V=(A1+A2+sqrt(A1A2))h/3,A1.A2为上下底面积,h为高

(以下为正棱台)

2. 侧面积 S=(p1+p2)l/2,p1.p2为上下底面周长,l为斜高

3. 全面积 T=S+A1+A2

圆柱:

1. 侧面积 S=2PIrh

2. 全面积 T=2PIr(h+r)

3. 体积 V=PIr^2h

圆锥:

1. 母线 l=sqrt(h^2+r^2)

2. 侧面积 S=PIrl

3. 全面积 T=PIr(l+r)

4. 体积 V=PIr^2h/3

圆台:

1. 母线 l=sqrt(h^2+(r1-r2)^2)

2. 侧面积 S=PI(r1+r2)l

3. 全面积 T=PIr1(l+r1)+PIr2(l+r2)

4. 体积 V=PI(r1^2+r2^2+r1r2)h/3

球:

1. 全面积 T=4PIr^2

2. 体积 V=4PIr^3/3

球台:

1. 侧面积 S=2PIrh

2. 全面积 T=PI(2rh+r1^2+r2^2)

3. 体积 V=PIh(3(r1^2+r2^2)+h^2)/6

球扇形:

1. 全面积 T=PIr(2h+r0),h为球冠高,r0为球冠底面半径

2. 体积 V=2PIr^2h/3

1.3 多边形

[cpp]  view plain  copy

1. #include <stdlib.h>  
2. #include <math.h>  
3. #define MAXN 1000  
4. #define offset 10000  
5. #define eps 1e-8  
6. #define zero(x) (((x)>0?(x):-(x))<eps)  
7. #define _sign(x) ((x)>eps?1:((x)<-eps?2:0))  
8. struct point{double x,y;};  
9. struct line{point a,b;};  
10.   
11.   
12. double xmult(point p1,point p2,point p0){  
13.     return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);  
14. }  
15.   
16.   
17. //判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线  
18. int is_convex(int n,point* p){  
19.     int i,s[3]={1,1,1};  
20.     for (i=0;i<n&&s[1]|s[2];i++)  
21.         s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0;  
22.     return s[1]|s[2];  
23. }  
24.   
25.   
26. //判定凸多边形,顶点按顺时针或逆时针给出,不允许相邻边共线  
27. int is_convex_v2(int n,point* p){  
28.     int i,s[3]={1,1,1};  
29.     for (i=0;i<n&&s[0]&&s[1]|s[2];i++)  
30.         s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0;  
31.     return s[0]&&s[1]|s[2];  
32. }  
33.   
34.   
35. //判点在凸多边形内或多边形边上,顶点按顺时针或逆时针给出  
36. int inside_convex(point q,int n,point* p){  
37.     int i,s[3]={1,1,1};  
38.     for (i=0;i<n&&s[1]|s[2];i++)  
39.         s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0;  
40.     return s[1]|s[2];  
41. }  
42.   
43.   
44. //判点在凸多边形内,顶点按顺时针或逆时针给出,在多边形边上返回0  
45. int inside_convex_v2(point q,int n,point* p){  
46.     int i,s[3]={1,1,1};  
47.     for (i=0;i<n&&s[0]&&s[1]|s[2];i++)  
48.         s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0;  
49.     return s[0]&&s[1]|s[2];  
50. }  
51.   
52.   
53. //判点在任意多边形内,顶点按顺时针或逆时针给出  
54. //on_edge表示点在多边形边上时的返回值,offset为多边形坐标上限  
55. int inside_polygon(point q,int n,point* p,int on_edge=1){  
56.     point q2;  
57.     int i=0,count;  
58.     while (i<n)  
59.         for (count=i=0,q2.x=rand()+offset,q2.y=rand()+offset;i<n;i++)  
60.             if (zero(xmult(q,p[i],p[(i+1)%n]))&&(p[i].x-q.x)*(p[(i+1)%n].x-q.x)<eps&&(p[i].y-q.y)*(p[(i+1)%n].y-q.y)<eps)  
61.                 return on_edge;  
62.             else if (zero(xmult(q,q2,p[i])))  
63.                 break;  
64.             else if (xmult(q,p[i],q2)*xmult(q,p[(i+1)%n],q2)<-eps&&xmult(p[i],q,p[(i+1)%n])*xmult(p[i],q2,p[(i+1)%n])<-eps)  
65.                 count++;  
66.     return count&1;  
67. }  
68.   
69.   
70. inline int opposite_side(point p1,point p2,point l1,point l2){  
71.     return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps;  
72. }  
73.   
74.   
75. inline int dot_online_in(point p,point l1,point l2){  
76.     return zero(xmult(p,l1,l2))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps;  
77. }  
78.   
79.   
80. //判线段在任意多边形内,顶点按顺时针或逆时针给出,与边界相交返回1  
81. int inside_polygon(point l1,point l2,int n,point* p){  
82.     point t[MAXN],tt;  
83.     int i,j,k=0;  
84.     if (!inside_polygon(l1,n,p)||!inside_polygon(l2,n,p))  
85.         return 0;  
86.     for (i=0;i<n;i++)  
87.         if (opposite_side(l1,l2,p[i],p[(i+1)%n])&&opposite_side(p[i],p[(i+1)%n],l1,l2))  
88.             return 0;  
89.         else if (dot_online_in(l1,p[i],p[(i+1)%n]))  
90.             t[k++]=l1;  
91.         else if (dot_online_in(l2,p[i],p[(i+1)%n]))  
92.             t[k++]=l2;  
93.         else if (dot_online_in(p[i],l1,l2))  
94.             t[k++]=p[i];  
95.     for (i=0;i<k;i++)  
96.         for (j=i+1;j<k;j++){  
97.             tt.x=(t[i].x+t[j].x)/2;  
98.             tt.y=(t[i].y+t[j].y)/2;  
99.             if (!inside_polygon(tt,n,p))  
100.                 return 0;             
101.         }  
102.     return 1;  
103. }  
104.   
105.   
106. point intersection(line u,line v){  
107.     point ret=u.a;  
108.     double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))  
109.             /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));  
110.     ret.x+=(u.b.x-u.a.x)*t;  
111.     ret.y+=(u.b.y-u.a.y)*t;  
112.     return ret;  
113. }  
114.   
115.   
116. point barycenter(point a,point b,point c){  
117.     line u,v;  
118.     u.a.x=(a.x+b.x)/2;  
119.     u.a.y=(a.y+b.y)/2;  
120.     u.b=c;  
121.     v.a.x=(a.x+c.x)/2;  
122.     v.a.y=(a.y+c.y)/2;  
123.     v.b=b;  
124.     return intersection(u,v);  
125. }  
126.   
127.   
128. //多边形重心  
129. point barycenter(int n,point* p){  
130.     point ret,t;  
131.     double t1=0,t2;  
132.     int i;  
133.     ret.x=ret.y=0;  
134.     for (i=1;i<n-1;i++)  
135.         if (fabs(t2=xmult(p[0],p[i],p[i+1]))>eps){  
136.             t=barycenter(p[0],p[i],p[i+1]);  
137.             ret.x+=t.x*t2;  
138.             ret.y+=t.y*t2;  
139.             t1+=t2;  
140.         }  
141.     if (fabs(t1)>eps)  
142.         ret.x/=t1,ret.y/=t1;  
143.     return ret;  
144. }  

1.4多边形切割

[cpp]  view plain  copy

1. //多边形切割  
2. //可用于半平面交  
3. #define MAXN 100  
4. #define eps 1e-8  
5. #define zero(x) (((x)>0?(x):-(x))<eps)  
6. struct point{double x,y;};  
7.   
8.   
9. double xmult(point p1,point p2,point p0){  
10.     return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);  
11. }  
12.   
13.   
14. int same_side(point p1,point p2,point l1,point l2){  
15.     return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps;  
16. }  
17.   
18.   
19. point intersection(point u1,point u2,point v1,point v2){  
20.     point ret=u1;  
21.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
22.             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
23.     ret.x+=(u2.x-u1.x)*t;  
24.     ret.y+=(u2.y-u1.y)*t;  
25.     return ret;  
26. }  
27.   
28.   
29. //将多边形沿l1,l2确定的直线切割在side侧切割,保证l1,l2,side不共线  
30. void polygon_cut(int& n,point* p,point l1,point l2,point side){  
31.     point pp[MAXN];  
32.     int m=0,i;  
33.     for (i=0;i<n;i++){  
34.         if (same_side(p[i],side,l1,l2))  
35.             pp[m++]=p[i];  
36.         if (!same_side(p[i],p[(i+1)%n],l1,l2)&&!(zero(xmult(p[i],l1,l2))&&zero(xmult(p[(i+1)%n],l1,l2))))  
37.             pp[m++]=intersection(p[i],p[(i+1)%n],l1,l2);  
38.     }  
39.     for (n=i=0;i<m;i++)  
40.         if (!i||!zero(pp[i].x-pp[i-1].x)||!zero(pp[i].y-pp[i-1].y))  
41.             p[n++]=pp[i];  
42.     if (zero(p[n-1].x-p[0].x)&&zero(p[n-1].y-p[0].y))  
43.         n--;  
44.     if (n<3)  
45.         n=0;  
46. }  

1.5 浮点函数

[cpp]  view plain  copy

1. //浮点几何函数库  
2. #include <math.h>  
3. #define eps 1e-8  
4. #define zero(x) (((x)>0?(x):-(x))<eps)  
5. struct point{double x,y;};  
6. struct line{point a,b;};  
7.   
8.   
9. //计算cross product (P1-P0)x(P2-P0)  
10. double xmult(point p1,point p2,point p0){  
11.     return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);  
12. }  
13. double xmult(double x1,double y1,double x2,double y2,double x0,double y0){  
14.     return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0);  
15. }  
16.   
17.   
18. //计算dot product (P1-P0).(P2-P0)  
19. double dmult(point p1,point p2,point p0){  
20.     return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y);  
21. }  
22. double dmult(double x1,double y1,double x2,double y2,double x0,double y0){  
23.     return (x1-x0)*(x2-x0)+(y1-y0)*(y2-y0);  
24. }  
25.   
26.   
27. //两点距离  
28. double distance(point p1,point p2){  
29.     return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));  
30. }  
31. double distance(double x1,double y1,double x2,double y2){  
32.     return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));  
33. }  
34.   
35.   
36. //判三点共线  
37. int dots_inline(point p1,point p2,point p3){  
38.     return zero(xmult(p1,p2,p3));  
39. }  
40. int dots_inline(double x1,double y1,double x2,double y2,double x3,double y3){  
41.     return zero(xmult(x1,y1,x2,y2,x3,y3));  
42. }  
43.   
44.   
45. //判点是否在线段上,包括端点  
46. int dot_online_in(point p,line l){  
47.     return zero(xmult(p,l.a,l.b))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps;  
48. }  
49. int dot_online_in(point p,point l1,point l2){  
50.     return zero(xmult(p,l1,l2))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps;  
51. }  
52. int dot_online_in(double x,double y,double x1,double y1,double x2,double y2){  
53.     return zero(xmult(x,y,x1,y1,x2,y2))&&(x1-x)*(x2-x)<eps&&(y1-y)*(y2-y)<eps;  
54. }  
55.   
56.   
57. //判点是否在线段上,不包括端点  
58. int dot_online_ex(point p,line l){  
59.     return dot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y))&&(!zero(p.x-l.b.x)||!zero(p.y-l.b.y));  
60. }  
61. int dot_online_ex(point p,point l1,point l2){  
62.     return dot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y))&&(!zero(p.x-l2.x)||!zero(p.y-l2.y));  
63. }  
64. int dot_online_ex(double x,double y,double x1,double y1,double x2,double y2){  
65.     return dot_online_in(x,y,x1,y1,x2,y2)&&(!zero(x-x1)||!zero(y-y1))&&(!zero(x-x2)||!zero(y-y2));  
66. }  
67.   
68.   
69. //判两点在线段同侧,点在线段上返回0  
70. int same_side(point p1,point p2,line l){  
71.     return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)>eps;  
72. }  
73. int same_side(point p1,point p2,point l1,point l2){  
74.     return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps;  
75. }  
76.   
77.   
78. //判两点在线段异侧,点在线段上返回0  
79. int opposite_side(point p1,point p2,line l){  
80.     return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps;  
81. }  
82. int opposite_side(point p1,point p2,point l1,point l2){  
83.     return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps;  
84. }  
85.   
86.   
87. //判两直线平行  
88. int parallel(line u,line v){  
89.     return zero((u.a.x-u.b.x)*(v.a.y-v.b.y)-(v.a.x-v.b.x)*(u.a.y-u.b.y));  
90. }  
91. int parallel(point u1,point u2,point v1,point v2){  
92.     return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y));  
93. }  
94.   
95.   
96. //判两直线垂直  
97. int perpendicular(line u,line v){  
98.     return zero((u.a.x-u.b.x)*(v.a.x-v.b.x)+(u.a.y-u.b.y)*(v.a.y-v.b.y));  
99. }  
100. int perpendicular(point u1,point u2,point v1,point v2){  
101.     return zero((u1.x-u2.x)*(v1.x-v2.x)+(u1.y-u2.y)*(v1.y-v2.y));  
102. }  
103.   
104.   
105. //判两线段相交,包括端点和部分重合  
106. int intersect_in(line u,line v){  
107.     if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))  
108.         return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);  
109.     return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);  
110. }  
111. int intersect_in(point u1,point u2,point v1,point v2){  
112.     if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2))  
113.         return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2);  
114.     return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2);  
115. }  
116.   
117.   
118. //判两线段相交,不包括端点和部分重合  
119. int intersect_ex(line u,line v){  
120.     return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);  
121. }  
122. int intersect_ex(point u1,point u2,point v1,point v2){  
123.     return opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2);  
124. }  
125.   
126.   
127. //计算两直线交点,注意事先判断直线是否平行!  
128. //线段交点请另外判线段相交(同时还是要判断是否平行!)  
129. point intersection(line u,line v){  
130.     point ret=u.a;  
131.     double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))  
132.             /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));  
133.     ret.x+=(u.b.x-u.a.x)*t;  
134.     ret.y+=(u.b.y-u.a.y)*t;  
135.     return ret;  
136. }  
137. point intersection(point u1,point u2,point v1,point v2){  
138.     point ret=u1;  
139.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
140.             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
141.     ret.x+=(u2.x-u1.x)*t;  
142.     ret.y+=(u2.y-u1.y)*t;  
143.     return ret;  
144. }  
145.   
146.   
147. //点到直线上的最近点  
148. point ptoline(point p,line l){  
149.     point t=p;  
150.     t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;  
151.     return intersection(p,t,l.a,l.b);  
152. }  
153. point ptoline(point p,point l1,point l2){  
154.     point t=p;  
155.     t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;  
156.     return intersection(p,t,l1,l2);  
157. }  
158.   
159.   
160. //点到直线距离  
161. double disptoline(point p,line l){  
162.     return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b);  
163. }  
164. double disptoline(point p,point l1,point l2){  
165.     return fabs(xmult(p,l1,l2))/distance(l1,l2);  
166. }  
167. double disptoline(double x,double y,double x1,double y1,double x2,double y2){  
168.     return fabs(xmult(x,y,x1,y1,x2,y2))/distance(x1,y1,x2,y2);  
169. }  
170.   
171.   
172. //点到线段上的最近点  
173. point ptoseg(point p,line l){  
174.     point t=p;  
175.     t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;  
176.     if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps)  
177.         return distance(p,l.a)<distance(p,l.b)?l.a:l.b;  
178.     return intersection(p,t,l.a,l.b);  
179. }  
180. point ptoseg(point p,point l1,point l2){  
181.     point t=p;  
182.     t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;  
183.     if (xmult(l1,t,p)*xmult(l2,t,p)>eps)  
184.         return distance(p,l1)<distance(p,l2)?l1:l2;  
185.     return intersection(p,t,l1,l2);  
186. }  
187.   
188.   
189. //点到线段距离  
190. double disptoseg(point p,line l){  
191.     point t=p;  
192.     t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;  
193.     if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps)  
194.         return distance(p,l.a)<distance(p,l.b)?distance(p,l.a):distance(p,l.b);  
195.     return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b);  
196. }  
197. double disptoseg(point p,point l1,point l2){  
198.     point t=p;  
199.     t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;  
200.     if (xmult(l1,t,p)*xmult(l2,t,p)>eps)  
201.         return distance(p,l1)<distance(p,l2)?distance(p,l1):distance(p,l2);  
202.     return fabs(xmult(p,l1,l2))/distance(l1,l2);  
203. }  
204.   
205.   
206. //矢量V以P为顶点逆时针旋转angle并放大scale倍  
207. point rotate(point v,point p,double angle,double scale){  
208.     point ret=p;  
209.     v.x-=p.x,v.y-=p.y;  
210.     p.x=scale*cos(angle);  
211.     p.y=scale*sin(angle);  
212.     ret.x+=v.x*p.x-v.y*p.y;  
213.     ret.y+=v.x*p.y+v.y*p.x;  
214.     return ret;  
215. }  
216.   
217.   
218. //p点关于直线L的对称点  
219. ponit symmetricalPointofLine(point p, line L)  
220. {  
221.     point p2;  
222.     double d;  
223.     d = L.a * L.a + L.b * L.b;  
224.     p2.x = (L.b * L.b * p.x - L.a * L.a * p.x -   
225.             2 * L.a * L.b * p.y - 2 * L.a * L.c) / d;  
226.     p2.y = (L.a * L.a * p.y - L.b * L.b * p.y -   
227.             2 * L.a * L.b * p.x - 2 * L.b * L.c) / d;  
228.     return p2;  
229. }  
230.   
231.   
232. //求两点的平分线  
233. line bisector(point& a, point& b) {  
234.     line ab, ans;  ab.set(a, b);  
235.     double midx = (a.x + b.x)/2.0,  midy = (a.y + b.y)/2.0;  
236.     ans.a = -ab.b, ans.b = -ab.a, ans.c = -ab.b * midx + ab.a * midy;  
237.     return ans;  
238. }  
239.   
240.   
241. // 已知入射线、镜面，求反射线。   
242. // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数;    
243. a2,b2,c2为入射光直线方程系数;    
244. a,b,c为反射光直线方程系数.   
245. // 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>.   
246. // 不要忘记结果中可能会有"negative zeros"   
247.   
248.   
249. void reflect(double a1,double b1,double c1,  
250. double a2,double b2,double c2,  
251. double &a,double &b,double &c)   
252. {   
253.     double n,m;   
254.     double tpb,tpa;   
255.     tpb=b1*b2+a1*a2;   
256.     tpa=a2*b1-a1*b2;   
257.     m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1);   
258.     n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1);   
259.     if(fabs(a1*b2-a2*b1)<1e-20)   
260.     {   
261.         a=a2;b=b2;c=c2;   
262.         return;   
263.     }   
264.     double xx,yy; //(xx,yy)是入射线与镜面的交点。   
265.     xx=(b1*c2-b2*c1)/(a1*b2-a2*b1);   
266.     yy=(a2*c1-a1*c2)/(a1*b2-a2*b1);   
267.     a=n;   
268.     b=-m;   
269.     c=m*yy-xx*n;   
270. }  

1.6 面积

[cpp]  view plain  copy

1. #include <math.h>  
2. struct point{double x,y;};  
3.   
4.   
5. //计算cross product (P1-P0)x(P2-P0)  
6. double xmult(point p1,point p2,point p0){  
7.     return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);  
8. }  
9. double xmult(double x1,double y1,double x2,double y2,double x0,double y0){  
10.     return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0);  
11. }  
12.   
13.   
14. //计算三角形面积,输入三顶点  
15. double area_triangle(point p1,point p2,point p3){  
16.     return fabs(xmult(p1,p2,p3))/2;  
17. }  
18. double area_triangle(double x1,double y1,double x2,double y2,double x3,double y3){  
19.     return fabs(xmult(x1,y1,x2,y2,x3,y3))/2;  
20. }  
21.   
22.   
23. //计算三角形面积,输入三边长  
24. double area_triangle(double a,double b,double c){  
25.     double s=(a+b+c)/2;  
26.     return sqrt(s*(s-a)*(s-b)*(s-c));  
27. }  
28.   
29.   
30. //计算多边形面积,顶点按顺时针或逆时针给出  
31. double area_polygon(int n,point* p){  
32.     double s1=0,s2=0;  
33.     int i;  
34.     for (i=0;i<n;i++)  
35.         s1+=p[(i+1)%n].y*p[i].x,s2+=p[(i+1)%n].y*p[(i+2)%n].x;  
36.     return fabs(s1-s2)/2;  
37. }  

1.7球面

[cpp]  view plain  copy

1. #include <math.h>  
2. const double pi=acos(-1);  
3.   
4.   
5. //计算圆心角lat表示纬度,-90<=w<=90,lng表示经度  
6. //返回两点所在大圆劣弧对应圆心角,0<=angle<=pi  
7. double angle(double lng1,double lat1,double lng2,double lat2){  
8.     double dlng=fabs(lng1-lng2)*pi/180;  
9.     while (dlng>=pi+pi)  
10.         dlng-=pi+pi;  
11.     if (dlng>pi)  
12.         dlng=pi+pi-dlng;  
13.     lat1*=pi/180,lat2*=pi/180;  
14.     return acos(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2));  
15. }  
16.   
17.   
18. //计算距离,r为球半径  
19. double line_dist(double r,double lng1,double lat1,double lng2,double lat2){  
20.     double dlng=fabs(lng1-lng2)*pi/180;  
21.     while (dlng>=pi+pi)  
22.         dlng-=pi+pi;  
23.     if (dlng>pi)  
24.         dlng=pi+pi-dlng;  
25.     lat1*=pi/180,lat2*=pi/180;  
26.     return r*sqrt(2-2*(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2)));  
27. }  
28.   
29.   
30. //计算球面距离,r为球半径  
31. inline double sphere_dist(double r,double lng1,double lat1,double lng2,double lat2){  
32.     return r*angle(lng1,lat1,lng2,lat2);  
33. }  
34.   
35.   
36. //球面反射  
37. //SGU110  
38. // http://acm.sgu.ru/problem.php?contest=0&problem=110  
39.   
40.   
41. #include <cstdio>  
42. #include <cmath>  
43.   
44.   
45. const int size = 555;  
46. const double eps = 1e-9;  
47.   
48.   
49. struct point {double x, y, z;} centre = {0, 0, 0};  
50. struct circle {point o; double r;} cir[size];  
51. struct ray {point s, dir;} l;  
52. int n;  
53.   
54.   
55. int dcmp (double x){return x < -eps ? -1 : x > eps;}  
56. double sqr (double x){return x*x;}  
57. double dot (point a, point b){return a.x * b.x + a.y * b.y + a.z * b.z;}  
58. double dis2 (point a, point b){return sqr(a.x-b.x) + sqr(a.y-b.y) + sqr(a.z-b.z);}  
59. double disToLine2 (point a, ray l){   /**** 点到直线L的距离的平方 **/  
60.     point tmp;  
61.     tmp.x =  l.dir.y * (a.z - l.s.z) - l.dir.z * (a.y - l.s.y);  
62.     tmp.y = -l.dir.x * (a.z - l.s.z) + l.dir.z * (a.x - l.s.x);  
63.     tmp.z =  l.dir.x * (a.y - l.s.y) - l.dir.y * (a.x - l.s.x);   
64.     return dis2 (tmp, centre) / dis2 (l.dir, centre);  
65. }  
66.   
67.   
68. /**　用解方程(点到圆心的距离为r)法求交点  (下面有向量法求交点, 两者取其一, 都OK)*/  
69. /* 是向量分量表示发的系数, 必须在射线上,故K必须为正, t是交点***/  
70. /* 
71. bool find (circle p, ray l, double &k, point &t)  
72. { 
73.     double x = l.s.x - p.o.x, y = l.s.y - p.o.y, z = l.s.z - p.o.z; 
74.     double a = sqr(l.dir.x) + sqr(l.dir.y) + sqr(l.dir.z); 
75.     double b = 2 * (x*l.dir.x + y*l.dir.y + z*l.dir.z); 
76.     double c = x*x + y*y + z*z - p.r*p.r; 
77.     double det = b*b - 4*a*c; 
78. //  printf ("a = %lf, b = %lf, c = %lf", a, b, c); 
79. //  printf ("det = %lf\n", det); 
80.     if (dcmp(det) == -1) return false; 
81.     k = (-b - sqrt(det)) / a / 2; 
82.     if (dcmp(k) != 1) return false; 
83.     t.x = l.s.x + k * l.dir.x; 
84.     t.y = l.s.y + k * l.dir.y; 
85.     t.z = l.s.z + k * l.dir.z; 
86.     return true; 
87. } 
88. */  
89.   
90.   
91. /**** 用向量法求交点  ***/  
92. bool find (circle p, ray l, double &k, point &t)  
93. {  
94.     double h2 = disToLine2 (p.o, l);  
95. //  printf ("h2 = %lf\n", h2);  
96.     if (dcmp(p.r*p.r - h2) < 0) return false;  
97.     point tmp;  
98.     tmp.x = p.o.x - l.s.x;  
99.     tmp.y = p.o.y - l.s.y;  
100.     tmp.z = p.o.z - l.s.z;  
101.     if (dcmp(dot(tmp, l.dir)) <= 0) return false;  
102.     k = sqrt(dis2(p.o, l.s) - h2) - sqrt(p.r*p.r - h2);  
103.     double k1 = k / sqrt(dis2(l.dir, centre));  
104.     t.x = l.s.x + k1 * l.dir.x;  
105.     t.y = l.s.y + k1 * l.dir.y;  
106.     t.z = l.s.z + k1 * l.dir.z;  
107.     return true;   
108. }  
109. /*计算新射线的起点和方向 */  
110. void newRay (ray &l, ray l1, point inter)  
111. {  
112.     double k = - 2 * dot(l.dir, l1.dir);  
113.     l.dir.x += l1.dir.x * k;  
114.     l.dir.y += l1.dir.y * k;  
115.     l.dir.z += l1.dir.z * k;  
116.     l.s = inter;  
117. }  
118. /* 返回的是最先相交的球的编号,均不相交,返回-1 */  
119. int update ()  
120. {  
121.     int sign = -1, i;  
122.     double k = 1e100, tmp;  
123.     point inter, t;  
124.     for (i = 1; i <= n; i++){ //找到最先相交的球  
125.         if (!find (cir[i], l, tmp, t)) continue;  
126.         if (dcmp (tmp - k) < 0) k = tmp, inter = t, sign = i;  
127.     }  
128.     //ray 变向  
129.     if (sign == -1) return sign;  
130.     ray l1;  
131.     l1.s = cir[sign].o;  
132.     l1.dir.x = (inter.x - l1.s.x) / cir[sign].r;  
133.     l1.dir.y = (inter.y - l1.s.y) / cir[sign].r;  
134.     l1.dir.z = (inter.z - l1.s.z) / cir[sign].r;  
135.     newRay (l, l1, inter);  
136.     return sign;  
137. }  
138. int main ()  
139. {  
140. //  freopen ("in", "r", stdin);  
141.     int i;  
142.     scanf ("%d", &n);  
143.     for (i = 1; i <= n; i++) //输入空间的球位置  
144.         scanf ("%lf%lf%lf%lf", &cir[i].o.x, &cir[i].o.y, &cir[i].o.z, &cir[i].r);  
145.     scanf ("%lf%lf%lf%lf%lf%lf", &l.s.x, &l.s.y, &l.s.z, &l.dir.x, &l.dir.y, &l.dir.z);  
146.     for (i = 0; i <= 10; i++){ //最多输出十次相交的球的编号  
147.         int sign = update ();  
148.         if (sign == -1) break;  
149.         if (i == 0) printf ("%d", sign);  
150.         else if (i < 10) printf (" %d", sign);  
151.         else printf (" etc.");  
152.     }  
153.     puts ("");  
154. }  

1.8三角形

[cpp]  view plain  copy

1. #include <math.h>  
2. struct point{double x,y;};  
3. struct line{point a,b;};  
4.   
5.   
6. double distance(point p1,point p2){  
7.     return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));  
8. }  
9.   
10.   
11. point intersection(line u,line v){  
12.     point ret=u.a;  
13.     double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))  
14.             /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));  
15.     ret.x+=(u.b.x-u.a.x)*t;  
16.     ret.y+=(u.b.y-u.a.y)*t;  
17.     return ret;  
18. }  
19.   
20.   
21. //外心  
22. point circumcenter(point a,point b,point c){  
23.     line u,v;  
24.     u.a.x=(a.x+b.x)/2;  
25.     u.a.y=(a.y+b.y)/2;  
26.     u.b.x=u.a.x-a.y+b.y;  
27.     u.b.y=u.a.y+a.x-b.x;  
28.     v.a.x=(a.x+c.x)/2;  
29.     v.a.y=(a.y+c.y)/2;  
30.     v.b.x=v.a.x-a.y+c.y;  
31.     v.b.y=v.a.y+a.x-c.x;  
32.     return intersection(u,v);  
33. }  
34.   
35.   
36. //内心  
37. point incenter(point a,point b,point c){  
38.     line u,v;  
39.     double m,n;  
40.     u.a=a;  
41.     m=atan2(b.y-a.y,b.x-a.x);  
42.     n=atan2(c.y-a.y,c.x-a.x);  
43.     u.b.x=u.a.x+cos((m+n)/2);  
44.     u.b.y=u.a.y+sin((m+n)/2);  
45.     v.a=b;  
46.     m=atan2(a.y-b.y,a.x-b.x);  
47.     n=atan2(c.y-b.y,c.x-b.x);  
48.     v.b.x=v.a.x+cos((m+n)/2);  
49.     v.b.y=v.a.y+sin((m+n)/2);  
50.     return intersection(u,v);  
51. }  
52.   
53.   
54. //垂心  
55. point perpencenter(point a,point b,point c){  
56.     line u,v;  
57.     u.a=c;  
58.     u.b.x=u.a.x-a.y+b.y;  
59.     u.b.y=u.a.y+a.x-b.x;  
60.     v.a=b;  
61.     v.b.x=v.a.x-a.y+c.y;  
62.     v.b.y=v.a.y+a.x-c.x;  
63.     return intersection(u,v);  
64. }  
65.   
66.   
67. //重心  
68. //到三角形三顶点距离的平方和最小的点  
69. //三角形内到三边距离之积最大的点  
70. point barycenter(point a,point b,point c){  
71.     line u,v;  
72.     u.a.x=(a.x+b.x)/2;  
73.     u.a.y=(a.y+b.y)/2;  
74.     u.b=c;  
75.     v.a.x=(a.x+c.x)/2;  
76.     v.a.y=(a.y+c.y)/2;  
77.     v.b=b;  
78.     return intersection(u,v);  
79. }  
80.   
81.   
82. //费马点  
83. //到三角形三顶点距离之和最小的点  
84. point fermentpoint(point a,point b,point c){  
85.     point u,v;  
86.     double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y);  
87.     int i,j,k;  
88.     u.x=(a.x+b.x+c.x)/3;  
89.     u.y=(a.y+b.y+c.y)/3;  
90.     while (step>1e-10)  
91.         for (k=0;k<10;step/=2,k++)  
92.             for (i=-1;i<=1;i++)  
93.                 for (j=-1;j<=1;j++){  
94.                     v.x=u.x+step*i;  
95.                     v.y=u.y+step*j;  
96.                     if (distance(u,a)+distance(u,b)+distance(u,c)>distance(v,a)+distance(v,b)+distance(v,c))  
97.                         u=v;  
98.                 }  
99.     return u;  
100. }  
101.   
102.   
103. //求曲率半径 三角形内最大可围成面积  
104. #include<iostream>  
105.  #include<cmath>  
106.  using namespace std;  
107.  const double pi=3.14159265358979;  
108.  int main()  
109.  {  
110.     double a,b,c,d,p,s,r,ans,R,x,l; int T=0;  
111.     while(cin>>a>>b>>c>>d&&a+b+c+d)  
112.      {  
113.         T++;  
114.         l=a+b+c;  
115.         p=l/2;  
116.         s=sqrt(p*(p-a)*(p-b)*(p-c));  
117.         R= s /p;  
118.         if (d >= l)  ans = s;  
119.         else if(2*pi*R>=d) ans=d*d/(4*pi);  
120.         else  
121.         {  
122.             r = (l-d)/((l/R)-(2*pi));  
123.             x = r*r*s/(R*R);  
124.             ans = s - x + pi * r * r;    
125.         }  
126.         printf("Case %d: %.2lf\n",T,ans);  
127.      }  
128.      return 0;  
129.  }  

1.9三维几何

[cpp]  view plain  copy

1. //三维几何函数库  
2. #include <math.h>  
3. #define eps 1e-8  
4. #define zero(x) (((x)>0?(x):-(x))<eps)  
5. struct point3{double x,y,z;};  
6. struct line3{point3 a,b;};  
7. struct plane3{point3 a,b,c;};  
8.   
9.   
10. //计算cross product U x V  
11. point3 xmult(point3 u,point3 v){  
12.     point3 ret;  
13.     ret.x=u.y*v.z-v.y*u.z;  
14.     ret.y=u.z*v.x-u.x*v.z;  
15.     ret.z=u.x*v.y-u.y*v.x;  
16.     return ret;  
17. }  
18.   
19.   
20. //计算dot product U . V  
21. double dmult(point3 u,point3 v){  
22.     return u.x*v.x+u.y*v.y+u.z*v.z;  
23. }  
24.   
25.   
26. //矢量差 U - V  
27. point3 subt(point3 u,point3 v){  
28.     point3 ret;  
29.     ret.x=u.x-v.x;  
30.     ret.y=u.y-v.y;  
31.     ret.z=u.z-v.z;  
32.     return ret;  
33. }  
34.   
35.   
36. //取平面法向量  
37. point3 pvec(plane3 s){  
38.     return xmult(subt(s.a,s.b),subt(s.b,s.c));  
39. }  
40. point3 pvec(point3 s1,point3 s2,point3 s3){  
41.     return xmult(subt(s1,s2),subt(s2,s3));  
42. }  
43.   
44.   
45. //两点距离,单参数取向量大小  
46. double distance(point3 p1,point3 p2){  
47.     return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));  
48. }  
49.   
50.   
51. //向量大小  
52. double vlen(point3 p){  
53.     return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);  
54. }  
55.   
56.   
57. //判三点共线  
58. int dots_inline(point3 p1,point3 p2,point3 p3){  
59.     return vlen(xmult(subt(p1,p2),subt(p2,p3)))<eps;  
60. }  
61.   
62.   
63. //判四点共面  
64. int dots_onplane(point3 a,point3 b,point3 c,point3 d){  
65.     return zero(dmult(pvec(a,b,c),subt(d,a)));  
66. }  
67.   
68.   
69. //判点是否在线段上,包括端点和共线  
70. int dot_online_in(point3 p,line3 l){  
71.     return zero(vlen(xmult(subt(p,l.a),subt(p,l.b))))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&  
72.         (l.a.y-p.y)*(l.b.y-p.y)<eps&&(l.a.z-p.z)*(l.b.z-p.z)<eps;  
73. }  
74. int dot_online_in(point3 p,point3 l1,point3 l2){  
75.     return zero(vlen(xmult(subt(p,l1),subt(p,l2))))&&(l1.x-p.x)*(l2.x-p.x)<eps&&  
76.         (l1.y-p.y)*(l2.y-p.y)<eps&&(l1.z-p.z)*(l2.z-p.z)<eps;  
77. }  
78.   
79.   
80. //判点是否在线段上,不包括端点  
81. int dot_online_ex(point3 p,line3 l){  
82.     return dot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y)||!zero(p.z-l.a.z))&&  
83.         (!zero(p.x-l.b.x)||!zero(p.y-l.b.y)||!zero(p.z-l.b.z));  
84. }  
85. int dot_online_ex(point3 p,point3 l1,point3 l2){  
86.     return dot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y)||!zero(p.z-l1.z))&&  
87.         (!zero(p.x-l2.x)||!zero(p.y-l2.y)||!zero(p.z-l2.z));  
88. }  
89.   
90.   
91. //判点是否在空间三角形上,包括边界,三点共线无意义  
92. int dot_inplane_in(point3 p,plane3 s){  
93.     return zero(vlen(xmult(subt(s.a,s.b),subt(s.a,s.c)))-vlen(xmult(subt(p,s.a),subt(p,s.b)))-  
94.         vlen(xmult(subt(p,s.b),subt(p,s.c)))-vlen(xmult(subt(p,s.c),subt(p,s.a))));  
95. }  
96. int dot_inplane_in(point3 p,point3 s1,point3 s2,point3 s3){  
97.     return zero(vlen(xmult(subt(s1,s2),subt(s1,s3)))-vlen(xmult(subt(p,s1),subt(p,s2)))-  
98.         vlen(xmult(subt(p,s2),subt(p,s3)))-vlen(xmult(subt(p,s3),subt(p,s1))));  
99. }  
100.   
101.   
102. //判点是否在空间三角形上,不包括边界,三点共线无意义  
103. int dot_inplane_ex(point3 p,plane3 s){  
104.     return dot_inplane_in(p,s)&&vlen(xmult(subt(p,s.a),subt(p,s.b)))>eps&&  
105.         vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps;  
106. }  
107. int dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3){  
108.     return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&&  
109.         vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps;  
110. }  
111.   
112.   
113. //判两点在线段同侧,点在线段上返回0,不共面无意义  
114. int same_side(point3 p1,point3 p2,line3 l){  
115.     return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))>eps;  
116. }  
117. int same_side(point3 p1,point3 p2,point3 l1,point3 l2){  
118.     return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))>eps;  
119. }  
120.   
121.   
122. //判两点在线段异侧,点在线段上返回0,不共面无意义  
123. int opposite_side(point3 p1,point3 p2,line3 l){  
124.     return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))<-eps;  
125. }  
126. int opposite_side(point3 p1,point3 p2,point3 l1,point3 l2){  
127.     return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))<-eps;  
128. }  
129.   
130.   
131. //判两点在平面同侧,点在平面上返回0  
132. int same_side(point3 p1,point3 p2,plane3 s){  
133.     return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps;  
134. }  
135. int same_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){  
136.     return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))>eps;  
137. }  
138.   
139.   
140. //判两点在平面异侧,点在平面上返回0  
141. int opposite_side(point3 p1,point3 p2,plane3 s){  
142.     return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps;  
143. }  
144. int opposite_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){  
145.     return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))<-eps;  
146. }  
147.   
148.   
149. //判两直线平行  
150. int parallel(line3 u,line3 v){  
151.     return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))<eps;  
152. }  
153. int parallel(point3 u1,point3 u2,point3 v1,point3 v2){  
154.     return vlen(xmult(subt(u1,u2),subt(v1,v2)))<eps;  
155. }  
156.   
157.   
158. //判两平面平行  
159. int parallel(plane3 u,plane3 v){  
160.     return vlen(xmult(pvec(u),pvec(v)))<eps;  
161. }  
162. int parallel(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){  
163.     return vlen(xmult(pvec(u1,u2,u3),pvec(v1,v2,v3)))<eps;  
164. }  
165.   
166.   
167. //判直线与平面平行  
168. int parallel(line3 l,plane3 s){  
169.     return zero(dmult(subt(l.a,l.b),pvec(s)));  
170. }  
171. int parallel(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
172.     return zero(dmult(subt(l1,l2),pvec(s1,s2,s3)));  
173. }  
174.   
175.   
176. //判两直线垂直  
177. int perpendicular(line3 u,line3 v){  
178.     return zero(dmult(subt(u.a,u.b),subt(v.a,v.b)));  
179. }  
180. int perpendicular(point3 u1,point3 u2,point3 v1,point3 v2){  
181.     return zero(dmult(subt(u1,u2),subt(v1,v2)));  
182. }  
183.   
184.   
185. //判两平面垂直  
186. int perpendicular(plane3 u,plane3 v){  
187.     return zero(dmult(pvec(u),pvec(v)));  
188. }  
189. int perpendicular(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){  
190.     return zero(dmult(pvec(u1,u2,u3),pvec(v1,v2,v3)));  
191. }  
192.   
193.   
194. //判直线与平面平行  
195. int perpendicular(line3 l,plane3 s){  
196.     return vlen(xmult(subt(l.a,l.b),pvec(s)))<eps;  
197. }  
198. int perpendicular(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
199.     return vlen(xmult(subt(l1,l2),pvec(s1,s2,s3)))<eps;  
200. }  
201.   
202.   
203. //判两线段相交,包括端点和部分重合  
204. int intersect_in(line3 u,line3 v){  
205.     if (!dots_onplane(u.a,u.b,v.a,v.b))  
206.         return 0;  
207.     if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))  
208.         return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);  
209.     return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);  
210. }  
211. int intersect_in(point3 u1,point3 u2,point3 v1,point3 v2){  
212.     if (!dots_onplane(u1,u2,v1,v2))  
213.         return 0;  
214.     if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2))  
215.         return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2);  
216.     return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2);  
217. }  
218.   
219.   
220. //判两线段相交,不包括端点和部分重合  
221. int intersect_ex(line3 u,line3 v){  
222.     return dots_onplane(u.a,u.b,v.a,v.b)&&opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);  
223. }  
224. int intersect_ex(point3 u1,point3 u2,point3 v1,point3 v2){  
225.     return dots_onplane(u1,u2,v1,v2)&&opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2);  
226. }  
227.   
228.   
229. //判线段与空间三角形相交,包括交于边界和(部分)包含  
230. int intersect_in(line3 l,plane3 s){  
231.     return !same_side(l.a,l.b,s)&&!same_side(s.a,s.b,l.a,l.b,s.c)&&  
232.         !same_side(s.b,s.c,l.a,l.b,s.a)&&!same_side(s.c,s.a,l.a,l.b,s.b);  
233. }  
234. int intersect_in(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
235.     return !same_side(l1,l2,s1,s2,s3)&&!same_side(s1,s2,l1,l2,s3)&&  
236.         !same_side(s2,s3,l1,l2,s1)&&!same_side(s3,s1,l1,l2,s2);  
237. }  
238.   
239.   
240. //判线段与空间三角形相交,不包括交于边界和(部分)包含  
241. int intersect_ex(line3 l,plane3 s){  
242.     return opposite_side(l.a,l.b,s)&&opposite_side(s.a,s.b,l.a,l.b,s.c)&&  
243.         opposite_side(s.b,s.c,l.a,l.b,s.a)&&opposite_side(s.c,s.a,l.a,l.b,s.b);  
244. }  
245. int intersect_ex(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
246.     return opposite_side(l1,l2,s1,s2,s3)&&opposite_side(s1,s2,l1,l2,s3)&&  
247.         opposite_side(s2,s3,l1,l2,s1)&&opposite_side(s3,s1,l1,l2,s2);  
248. }  
249.   
250.   
251. //计算两直线交点,注意事先判断直线是否共面和平行!  
252. //线段交点请另外判线段相交(同时还是要判断是否平行!)  
253. point3 intersection(line3 u,line3 v){  
254.     point3 ret=u.a;  
255.     double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))  
256.             /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));  
257.     ret.x+=(u.b.x-u.a.x)*t;  
258.     ret.y+=(u.b.y-u.a.y)*t;  
259.     ret.z+=(u.b.z-u.a.z)*t;  
260.     return ret;  
261. }  
262. point3 intersection(point3 u1,point3 u2,point3 v1,point3 v2){  
263.     point3 ret=u1;  
264.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
265.             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
266.     ret.x+=(u2.x-u1.x)*t;  
267.     ret.y+=(u2.y-u1.y)*t;  
268.     ret.z+=(u2.z-u1.z)*t;  
269.     return ret;  
270. }  
271.   
272.   
273. //计算直线与平面交点,注意事先判断是否平行,并保证三点不共线!  
274. //线段和空间三角形交点请另外判断  
275. point3 intersection(line3 l,plane3 s){  
276.     point3 ret=pvec(s);  
277.     double t=(ret.x*(s.a.x-l.a.x)+ret.y*(s.a.y-l.a.y)+ret.z*(s.a.z-l.a.z))/  
278.         (ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z));  
279.     ret.x=l.a.x+(l.b.x-l.a.x)*t;  
280.     ret.y=l.a.y+(l.b.y-l.a.y)*t;  
281.     ret.z=l.a.z+(l.b.z-l.a.z)*t;  
282.     return ret;  
283. }  
284. point3 intersection(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
285.     point3 ret=pvec(s1,s2,s3);  
286.     double t=(ret.x*(s1.x-l1.x)+ret.y*(s1.y-l1.y)+ret.z*(s1.z-l1.z))/  
287.         (ret.x*(l2.x-l1.x)+ret.y*(l2.y-l1.y)+ret.z*(l2.z-l1.z));  
288.     ret.x=l1.x+(l2.x-l1.x)*t;  
289.     ret.y=l1.y+(l2.y-l1.y)*t;  
290.     ret.z=l1.z+(l2.z-l1.z)*t;  
291.     return ret;  
292. }  
293.   
294.   
295. //计算两平面交线,注意事先判断是否平行,并保证三点不共线!  
296. line3 intersection(plane3 u,plane3 v){  
297.     line3 ret;  
298.     ret.a=parallel(v.a,v.b,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.a,v.b,u.a,u.b,u.c);  
299.     ret.b=parallel(v.c,v.a,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.c,v.a,u.a,u.b,u.c);  
300.     return ret;  
301. }  
302. line3 intersection(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){  
303.     line3 ret;  
304.     ret.a=parallel(v1,v2,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v1,v2,u1,u2,u3);  
305.     ret.b=parallel(v3,v1,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v3,v1,u1,u2,u3);  
306.     return ret;  
307. }  
308.   
309.   
310. //点到直线距离  
311. double ptoline(point3 p,line3 l){  
312.     return vlen(xmult(subt(p,l.a),subt(l.b,l.a)))/distance(l.a,l.b);  
313. }  
314. double ptoline(point3 p,point3 l1,point3 l2){  
315.     return vlen(xmult(subt(p,l1),subt(l2,l1)))/distance(l1,l2);  
316. }  
317.   
318.   
319. //点到平面距离  
320. double ptoplane(point3 p,plane3 s){  
321.     return fabs(dmult(pvec(s),subt(p,s.a)))/vlen(pvec(s));  
322. }  
323. double ptoplane(point3 p,point3 s1,point3 s2,point3 s3){  
324.     return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));  
325. }  
326.   
327.   
328. //直线到直线距离  
329. double linetoline(line3 u,line3 v){  
330.     point3 n=xmult(subt(u.a,u.b),subt(v.a,v.b));  
331.     return fabs(dmult(subt(u.a,v.a),n))/vlen(n);  
332. }  
333. double linetoline(point3 u1,point3 u2,point3 v1,point3 v2){  
334.     point3 n=xmult(subt(u1,u2),subt(v1,v2));  
335.     return fabs(dmult(subt(u1,v1),n))/vlen(n);  
336. }  
337.   
338.   
339. //两直线夹角cos值  
340. double angle_cos(line3 u,line3 v){  
341.     return dmult(subt(u.a,u.b),subt(v.a,v.b))/vlen(subt(u.a,u.b))/vlen(subt(v.a,v.b));  
342. }  
343. double angle_cos(point3 u1,point3 u2,point3 v1,point3 v2){  
344.     return dmult(subt(u1,u2),subt(v1,v2))/vlen(subt(u1,u2))/vlen(subt(v1,v2));  
345. }  
346.   
347.   
348. //两平面夹角cos值  
349. double angle_cos(plane3 u,plane3 v){  
350.     return dmult(pvec(u),pvec(v))/vlen(pvec(u))/vlen(pvec(v));  
351. }  
352. double angle_cos(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){  
353.     return dmult(pvec(u1,u2,u3),pvec(v1,v2,v3))/vlen(pvec(u1,u2,u3))/vlen(pvec(v1,v2,v3));  
354. }  
355.   
356.   
357. //直线平面夹角sin值  
358. double angle_sin(line3 l,plane3 s){  
359.     return dmult(subt(l.a,l.b),pvec(s))/vlen(subt(l.a,l.b))/vlen(pvec(s));  
360. }  
361. double angle_sin(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){  
362.     return dmult(subt(l1,l2),pvec(s1,s2,s3))/vlen(subt(l1,l2))/vlen(pvec(s1,s2,s3));  
363. }  

1.10 凸包

[cpp]  view plain  copy

1. //水平序  
2. #define maxn 100005  
3.   
4.   
5. struct point  
6. {double x,y;}p[maxn],s[maxn];  
7. bool operator < (point a,point b)  
8. {return a.x<b.x || a.x==b.x&&a.y<b.y;}  
9.   
10.   
11. int n,f;  
12.   
13.   
14. double cp(point a,point b,point c)  
15. {return (c.y-a.y)*(b.x-a.x)-(b.y-a.y)*(c.x-a.x);}  
16.   
17.   
18. void Convex(point *p,int &n)  
19. {  
20.     sort(p,p+n);  
21.     int i,j,r,top,m;  
22.     s[0] = p[0];s[1] = p[1];top = 1;  
23.     for(i=2;i<n;i++)  
24.     {  
25.         while( top>0 && cp(p[i],s[top],s[top-1])>=0 ) top--;  
26.         top++;s[top] = p[i];  
27.     }  
28.     m = top;  
29.     top++;s[top] = p[n-2];  
30.     for(i=n-3;i>=0;i--)  
31.     {  
32.         while( top>m && cp(p[i],s[top],s[top-1])>=0 ) top--;  
33.         top++;s[top] = p[i];  
34.     }  
35.     top--;  
36.     n = top+1;  
37. }  
38.   
39.   
40. 极角序  
41. #include <stdio.h>  
42. #include <string.h>  
43. #include <algorithm>  
44. #include <math.h>  
45. using namespace std;  
46. #define maxn 100005  
47. int N;  
48. struct A  
49. {  
50.     int x,y;  
51.     int v,l;  
52. }P[maxn];  
53. int xmult(int x1,int y1,int x2,int y2,int x3,int y3)  
54. {  
55.     return (y2-y1)*(x3-x1)-(y3-y1)*(x2-x1);  
56. }  
57. void swap(A &a,A &b)  
58. {  
59.     A t = a;a = b,b = t;  
60. }  
61. bool operator < (A a,A b)  
62. {  
63.     int k = xmult(P[0].x,P[0].y,a.x,a.y,b.x,b.y);  
64.     if( k<0 )  
65.         return 1;  
66.     else if( k==0 )  
67.     {  
68.         if( abs(P[0].x-a.x)<abs(P[0].x-b.x) )  
69.             return 1;  
70.         if( abs(P[0].y-a.y)<abs(P[0].y-b.y) )  
71.             return 1;  
72.     }  
73.     return 0;  
74. }  
75. void Grem_scan(int n)  
76. {  
77.     int i,j,k,l;  
78.     k = 0x7fffffff;  
79.     for(i=0;i<n;i++)  
80.         if( P[i].x<k || P[i].x==k && P[i].y<P[l].y )  
81.         k = P[i].x,l = i;  
82.     swap(P[l],P[0]);  
83.     sort(P+1,P+n);  
84.       
85.     l = 3;  
86.     for(i=3;i<n;i++)  
87.     {  
88.         while( xmult(P[l-2].x,P[l-2].y,P[l-1].x,P[l-1].y,P[i].x,P[i].y)>0 )  
89.             l--;  
90.         P[l++] = P[i];  
91.     }  
92. }  
93. main()  
94. {  
95.     int i,j,k,l;  
96.     N = 0;  
97.     while( scanf("%d%d",&P[N].x,&P[N].y)!=EOF )  
98.         N++;  
99.     Grem_scan(N);  
100.     for(i=0;i<N;i++)  
101.         if( P[i].x==0 && P[i].y==0 )  
102.         break;  
103.     k = i++;  
104.     printf("(0,0)\n");  
105.     while( i!=k )  
106.         printf("(%d,%d)\n",P[i].x,P[i].y),i = (i+1)%N;  
107. }  
108.   
109.   
110. //卷包裹法  
111. #include <stdio.h>  
112. #include <string.h>  
113. #include <algorithm>  
114. using namespace std;  
115. #define maxn 55  
116. struct A  
117. {  
118.     int x,y;  
119. }P[maxn];  
120. int T,N;  
121. bool B[maxn];  
122. int as[maxn],L;  
123. int xmult(A a,A b,A c)  
124. {  
125.     return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);  
126. }  
127. main()  
128. {  
129.     int i,j,k,l;  
130.     scanf("%d",&T);  
131.     while( T-- )  
132.     {  
133.         scanf("%d",&N);  
134.         k = 0x7ffffff;  
135.         for(i=0;i<N;i++)  
136.         {  
137.             scanf("%d%d%d",&j,&P[i].x,&P[i].y);  
138.             if( P[i].y<k )  
139.                 k = P[i].y,l = i;  
140.         }  
141.         memset(B,0,sizeof(B));  
142.         B[l] = 1;  
143.         as[0] = l;  
144.         L = 1;  
145.         while( 1 )  
146.         {  
147.             A a,b;  
148.             if( L==1 )  
149.                 a.x = 0,a.y = P[as[0]].y;  
150.             else  
151.                 a = P[as[L-2]];  
152.             b = P[as[L-1]];  
153.   
154.   
155.             k = -1;  
156.             for(i=0;i<N;i++)  
157.             {  
158.                 if( B[i] )  
159.                     continue;  
160.                 if( xmult(a,b,P[i])<0 )  
161.                     continue;  
162.                 if( k==-1 || xmult(P[as[L-1]],P[k],P[i])<0 || xmult(P[as[L-1]],P[k],P[i])==0 && P[i].y<P[k].y )  
163.                     k = i;  
164.             }  
165.             if( k==-1 )  
166.                 break;  
167.             B[k] = 1;  
168.             as[L++] = k;  
169.         }  
170.         printf("%d ",L);  
171.         for(i=0;i<L;i++)  
172.             printf("%d ",as[i]+1);  
173.         printf("\n");  
174.     }  
175. }  
176. }  

1.11 网格

[cpp]  view plain  copy

1. #define abs(x) ((x)>0?(x):-(x))  
2. struct point{int x,y;};  
3.   
4.   
5. int gcd(int a,int b){return b?gcd(b,a%b):a;}  
6.   
7.   
8. //多边形上的网格点个数  
9. int grid_onedge(int n,point* p){  
10.     int i,ret=0;  
11.     for (i=0;i<n;i++)  
12.         ret+=gcd(abs(p[i].x-p[(i+1)%n].x),abs(p[i].y-p[(i+1)%n].y));  
13.     return ret;  
14. }  
15.   
16.   
17. //多边形内的网格点个数  
18. int grid_inside(int n,point* p){  
19.     int i,ret=0;  
20.     for (i=0;i<n;i++)  
21.         ret+=p[(i+1)%n].y*(p[i].x-p[(i+2)%n].x);  
22.     return (abs(ret)-grid_onedge(n,p))/2+1;  
23. }  

1.12 圆

[cpp]  view plain  copy

1. #include <math.h>  
2. #define eps 1e-8  
3. struct point{double x,y;};  
4.   
5.   
6. double xmult(point p1,point p2,point p0){  
7.     return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);  
8. }  
9.   
10.   
11. double distance(point p1,point p2){  
12.     return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));  
13. }  
14.   
15.   
16. double disptoline(point p,point l1,point l2){  
17.     return fabs(xmult(p,l1,l2))/distance(l1,l2);  
18. }  
19.   
20.   
21. point intersection(point u1,point u2,point v1,point v2){  
22.     point ret=u1;  
23.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
24.             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
25.     ret.x+=(u2.x-u1.x)*t;  
26.     ret.y+=(u2.y-u1.y)*t;  
27.     return ret;  
28. }  
29.   
30.   
31. //判直线和圆相交,包括相切  
32. int intersect_line_circle(point c,double r,point l1,point l2){  
33.     return disptoline(c,l1,l2)<r+eps;  
34. }  
35.   
36.   
37. //判线段和圆相交,包括端点和相切  
38. int intersect_seg_circle(point c,double r,point l1,point l2){  
39.     double t1=distance(c,l1)-r,t2=distance(c,l2)-r;  
40.     point t=c;  
41.     if (t1<eps||t2<eps)  
42.         return t1>-eps||t2>-eps;  
43.     t.x+=l1.y-l2.y;  
44.     t.y+=l2.x-l1.x;  
45.     return xmult(l1,c,t)*xmult(l2,c,t)<eps&&disptoline(c,l1,l2)-r<eps;  
46. }  
47.   
48.   
49. //判圆和圆相交,包括相切  
50. int intersect_circle_circle(point c1,double r1,point c2,double r2){  
51.     return distance(c1,c2)<r1+r2+eps&&distance(c1,c2)>fabs(r1-r2)-eps;  
52. }  
53.   
54.   
55. //计算圆上到点p最近点,如p与圆心重合,返回p本身  
56. point dot_to_circle(point c,double r,point p){  
57.     point u,v;  
58.     if (distance(p,c)<eps)  
59.         return p;  
60.     u.x=c.x+r*fabs(c.x-p.x)/distance(c,p);  
61.     u.y=c.y+r*fabs(c.y-p.y)/distance(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1);  
62.     v.x=c.x-r*fabs(c.x-p.x)/distance(c,p);  
63.     v.y=c.y-r*fabs(c.y-p.y)/distance(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1);  
64.     return distance(u,p)<distance(v,p)?u:v;  
65. }  
66.   
67.   
68. //计算直线与圆的交点,保证直线与圆有交点  
69. //计算线段与圆的交点可用这个函数后判点是否在线段上  
70. void intersection_line_circle(point c,double r,point l1,point l2,point& p1,point& p2){  
71.     point p=c;  
72.     double t;  
73.     p.x+=l1.y-l2.y;  
74.     p.y+=l2.x-l1.x;  
75.     p=intersection(p,c,l1,l2);  
76.     t=sqrt(r*r-distance(p,c)*distance(p,c))/distance(l1,l2);  
77.     p1.x=p.x+(l2.x-l1.x)*t;  
78.     p1.y=p.y+(l2.y-l1.y)*t;  
79.     p2.x=p.x-(l2.x-l1.x)*t;  
80.     p2.y=p.y-(l2.y-l1.y)*t;  
81. }  
82.   
83.   
84. //计算圆与圆的交点,保证圆与圆有交点,圆心不重合  
85. void intersection_circle_circle(point c1,double r1,point c2,double r2,point& p1,point& p2){  
86.     point u,v;  
87.     double t;  
88.     t=(1+(r1*r1-r2*r2)/distance(c1,c2)/distance(c1,c2))/2;  
89.     u.x=c1.x+(c2.x-c1.x)*t;  
90.     u.y=c1.y+(c2.y-c1.y)*t;  
91.     v.x=u.x+c1.y-c2.y;  
92.     v.y=u.y-c1.x+c2.x;  
93.     intersection_line_circle(c1,r1,u,v,p1,p2);  
94. }  
95.   
96.   
97. //将向量p逆时针旋转angle角度  
98. Point Rotate(Point p,double angle) {  
99.     Point res;  
100.     res.x=p.x*cos(angle)-p.y*sin(angle);  
101.     res.y=p.x*sin(angle)+p.y*cos(angle);  
102.     return res;  
103. }  
104. //求圆外一点对圆(o,r)的两个切点result1和result2  
105. void TangentPoint_PC(Point poi,Point o,double r,Point &result1,Point &result2) {  
106.     double line=sqrt((poi.x-o.x)*(poi.x-o.x)+(poi.y-o.y)*(poi.y-o.y));  
107.     double angle=acos(r/line);  
108.     Point unitvector,lin;  
109.     lin.x=poi.x-o.x;  
110.     lin.y=poi.y-o.y;  
111.     unitvector.x=lin.x/sqrt(lin.x*lin.x+lin.y*lin.y)*r;  
112.     unitvector.y=lin.y/sqrt(lin.x*lin.x+lin.y*lin.y)*r;  
113.     result1=Rotate(unitvector,-angle);  
114.     result2=Rotate(unitvector,angle);  
115.     result1.x+=o.x;  
116.     result1.y+=o.y;  
117.     result2.x+=o.x;  
118.     result2.y+=o.y;  
119.     return;  
120. }  

1.13 矢量运算求几何模板

[cpp]  view plain  copy

1. #include <iostream>  
2. #include <cmath>   
3. #include <vector>   
4. #include <algorithm>   
5. #define MAX_N 100  
6. using namespace std;   
7.   
8.   
9. ///  
10. //常量区  
11. const double INF        = 1e10;     // 无穷大   
12. const double EPS        = 1e-15;    // 计算精度   
13. const int LEFT          = 0;        // 点在直线左边   
14. const int RIGHT         = 1;        // 点在直线右边   
15. const int ONLINE        = 2;        // 点在直线上   
16. const int CROSS         = 0;        // 两直线相交   
17. const int COLINE        = 1;        // 两直线共线   
18. const int PARALLEL      = 2;        // 两直线平行   
19. const int NOTCOPLANAR   = 3;        // 两直线不共面   
20. const int INSIDE        = 1;        // 点在图形内部   
21. const int OUTSIDE       = 2;        // 点在图形外部   
22. const int BORDER        = 3;        // 点在图形边界   
23. const int BAOHAN        = 1;        // 大圆包含小圆  
24. const int NEIQIE        = 2;        // 内切  
25. const int XIANJIAO      = 3;        // 相交  
26. const int WAIQIE        = 4;        // 外切  
27. const int XIANLI        = 5;        // 相离  
28. const double pi        = acos(-1.0)  //圆周率  
29. ///   
30.   
31.   
32. ///  
33. //类型定义区  
34. struct Point {              // 二维点或矢量   
35.     double x, y;   
36.     double angle, dis;   
37.     Point() {}   
38.     Point(double x0, double y0): x(x0), y(y0) {}   
39. };   
40. struct Point3D {            //三维点或矢量   
41.     double x, y, z;   
42.     Point3D() {}   
43.     Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {}   
44. };   
45. struct Line {               // 二维的直线或线段   
46.     Point p1, p2;   
47.     Line() {}   
48.     Line(Point p10, Point p20): p1(p10), p2(p20) {}   
49. };   
50. struct Line3D {             // 三维的直线或线段   
51.     Point3D p1, p2;   
52.     Line3D() {}   
53.     Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {}   
54. };   
55. struct Rect {              // 用长宽表示矩形的方法 w, h分别表示宽度和高度   
56.     double w, h;   
57.  Rect() {}  
58.  Rect(double _w,double _h) : w(_w),h(_h) {}  
59. };   
60. struct Rect_2 {             // 表示矩形，左下角坐标是(xl, yl)，右上角坐标是(xh, yh)   
61.     double xl, yl, xh, yh;   
62.  Rect_2() {}  
63.  Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {}  
64. };   
65. struct Circle {            //圆  
66.  Point c;  
67.  double r;  
68.  Circle() {}  
69.  Circle(Point _c,double _r) :c(_c),r(_r) {}  
70. };  
71. typedef vector<Point> Polygon;      // 二维多边形       
72. typedef vector<Point> Points;       // 二维点集   
73. typedef vector<Point3D> Points3D;   // 三维点集   
74. ///   
75.   
76.   
77. ///  
78. //基本函数区  
79. inline double max(double x,double y)   
80. {   
81.     return x > y ? x : y;   
82. }   
83. inline double min(double x, double y)   
84. {   
85.     return x > y ? y : x;   
86. }   
87. inline bool ZERO(double x)              // x == 0   
88. {   
89.     return (fabs(x) < EPS);   
90. }   
91. inline bool ZERO(Point p)               // p == 0   
92. {   
93.     return (ZERO(p.x) && ZERO(p.y));   
94. }   
95. inline bool ZERO(Point3D p)              // p == 0   
96. {   
97.     return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z));   
98. }   
99. inline bool EQ(double x, double y)      // eqaul, x == y   
100. {   
101.     return (fabs(x - y) < EPS);   
102. }   
103. inline bool NEQ(double x, double y)     // not equal, x != y   
104. {   
105.     return (fabs(x - y) >= EPS);   
106. }   
107. inline bool LT(double x, double y)     // less than, x < y   
108. {   
109.     return ( NEQ(x, y) && (x < y) );   
110. }   
111. inline bool GT(double x, double y)     // greater than, x > y   
112. {   
113.     return ( NEQ(x, y) && (x > y) );   
114. }   
115. inline bool LEQ(double x, double y)     // less equal, x <= y   
116. {   
117.     return ( EQ(x, y) || (x < y) );   
118. }   
119. inline bool GEQ(double x, double y)     // greater equal, x >= y   
120. {   
121.     return ( EQ(x, y) || (x > y) );   
122. }   
123. // 注意！！！   
124. // 如果是一个很小的负的浮点数   
125. // 保留有效位数输出的时候会出现-0.000这样的形式，   
126. // 前面多了一个负号   
127. // 这就会导致错误！！！！！！   
128. // 因此在输出浮点数之前，一定要调用次函数进行修正！   
129. inline double FIX(double x)   
130. {   
131.     return (fabs(x) < EPS) ? 0 : x;   
132. }   
133. //   
134.   
135.   
136. /  
137. //二维矢量运算   
138. bool operator==(Point p1, Point p2)    
139. {   
140.     return ( EQ(p1.x, p2.x) &&  EQ(p1.y, p2.y) );   
141. }   
142. bool operator!=(Point p1, Point p2)    
143. {   
144.     return ( NEQ(p1.x, p2.x) ||  NEQ(p1.y, p2.y) );   
145. }   
146. bool operator<(Point p1, Point p2)   
147. {   
148.     if (NEQ(p1.x, p2.x)) {   
149.         return (p1.x < p2.x);   
150.     } else {   
151.         return (p1.y < p2.y);   
152.     }   
153. }   
154. Point operator+(Point p1, Point p2)    
155. {   
156.     return Point(p1.x + p2.x, p1.y + p2.y);   
157. }   
158. Point operator-(Point p1, Point p2)    
159. {   
160.     return Point(p1.x - p2.x, p1.y - p2.y);   
161. }   
162. double operator*(Point p1, Point p2) // 计算叉乘 p1 × p2   
163. {   
164.     return (p1.x * p2.y - p2.x * p1.y);   
165. }   
166. double operator&(Point p1, Point p2) { // 计算点积 p1·p2   
167.     return (p1.x * p2.x + p1.y * p2.y);   
168. }   
169. double Norm(Point p) // 计算矢量p的模   
170. {   
171.     return sqrt(p.x * p.x + p.y * p.y);   
172. }   
173. // 把矢量p旋转角度angle (弧度表示)   
174. // angle > 0表示逆时针旋转   
175. // angle < 0表示顺时针旋转   
176. Point Rotate(Point p, double angle)   
177. {   
178.     Point result;   
179.     result.x = p.x * cos(angle) - p.y * sin(angle);   
180.     result.y = p.x * sin(angle) + p.y * cos(angle);   
181.     return result;   
182. }   
183. //   
184.   
185.   
186. //   
187. //三维矢量运算   
188. bool operator==(Point3D p1, Point3D p2)    
189. {   
190.     return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) && EQ(p1.z, p2.z) );   
191. }   
192. bool operator<(Point3D p1, Point3D p2)   
193. {   
194.     if (NEQ(p1.x, p2.x)) {   
195.         return (p1.x < p2.x);   
196.     } else if (NEQ(p1.y, p2.y)) {   
197.         return (p1.y < p2.y);   
198.     } else {   
199.         return (p1.z < p2.z);   
200.     }   
201. }   
202. Point3D operator+(Point3D p1, Point3D p2)    
203. {   
204.     return Point3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z);   
205. }   
206. Point3D operator-(Point3D p1, Point3D p2)    
207. {   
208.     return Point3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z);   
209. }   
210. Point3D operator*(Point3D p1, Point3D p2) // 计算叉乘 p1 x p2   
211. {   
212.     return Point3D(p1.y * p2.z - p1.z * p2.y,   
213.         p1.z * p2.x - p1.x * p2.z,   
214.         p1.x * p2.y - p1.y * p2.x );           
215. }   
216. double operator&(Point3D p1, Point3D p2) { // 计算点积 p1·p2   
217.     return (p1.x * p2.x + p1.y * p2.y + p1.z * p2.z);   
218. }   
219. double Norm(Point3D p) // 计算矢量p的模   
220. {   
221.     return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);   
222. }   
223.   
224.   
225. //   
226.   
227.   
228. /  
229. //点.线段.直线问题  
230. //  
231. double Distance(Point p1, Point p2) //2点间的距离  
232. {  
233.  return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));  
234. }  
235. double Distance(Point3D p1, Point3D p2) //2点间的距离,三维  
236. {  
237.  return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));  
238. }  
239. double Distance(Point p, Line L) // 求二维平面上点到直线的距离   
240. {   
241.     return ( fabs((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) );   
242. }   
243. double Distance(Point3D p, Line3D L)// 求三维空间中点到直线的距离   
244. {   
245.     return ( Norm((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) );   
246. }   
247. bool OnLine(Point p, Line L) // 判断二维平面上点p是否在直线L上   
248. {   
249.     return ZERO( (p - L.p1) * (L.p2 - L.p1) );   
250. }   
251. bool OnLine(Point3D p, Line3D L) // 判断三维空间中点p是否在直线L上   
252. {   
253.     return ZERO( (p - L.p1) * (L.p2 - L.p1) );   
254. }   
255. int Relation(Point p, Line L) // 计算点p与直线L的相对关系 ,返回ONLINE,LEFT,RIGHT  
256. {   
257.     double res = (L.p2 - L.p1) * (p - L.p1);   
258.     if (EQ(res, 0)) {   
259.         return ONLINE;   
260.     } else if (res > 0) {   
261.         return LEFT;   
262.     } else {   
263.         return RIGHT;   
264.     }   
265. }   
266. bool SameSide(Point p1, Point p2, Line L) // 判断点p1, p2是否在直线L的同侧   
267. {   
268.     double m1 = (p1 - L.p1) * (L.p2 - L.p1);   
269.     double m2 = (p2 - L.p1) * (L.p2 - L.p1);   
270.     return GT(m1 * m2, 0);   
271. }   
272. bool OnLineSeg(Point p, Line L) // 判断二维平面上点p是否在线段l上   
273. {   
274.     return ( ZERO( (L.p1 - p) * (L.p2 - p) ) &&   
275.         LEQ((p.x - L.p1.x)*(p.x - L.p2.x), 0) &&   
276.         LEQ((p.y - L.p1.y)*(p.y - L.p2.y), 0) );   
277. }   
278. bool OnLineSeg(Point3D p, Line3D L) // 判断三维空间中点p是否在线段l上   
279. {   
280.     return ( ZERO((L.p1 - p) * (L.p2 - p)) &&   
281.         EQ( Norm(p - L.p1) + Norm(p - L.p2), Norm(L.p2 - L.p1)) );   
282. }   
283. Point SymPoint(Point p, Line L) // 求二维平面上点p关于直线L的对称点   
284. {   
285.     Point result;   
286.     double a = L.p2.x - L.p1.x;   
287.     double b = L.p2.y - L.p1.y;   
288.     double t = ( (p.x - L.p1.x) * a + (p.y - L.p1.y) * b ) / (a*a + b*b);   
289.     result.x = 2 * L.p1.x + 2 * a * t - p.x;   
290.     result.y = 2 * L.p1.y + 2 * b * t - p.y;   
291.     return result;   
292. }   
293. bool Coplanar(Points3D points) // 判断一个点集中的点是否全部共面   
294. {   
295.     int i;   
296.     Point3D p;   
297.   
298.   
299.     if (points.size() < 4) return true;   
300.     p = (points[2] - points[0]) * (points[1] - points[0]);   
301.     for (i = 3; i < points.size(); i++) {   
302.         if (! ZERO(p & points[i]) ) return false;   
303.     }   
304.     return true;   
305. }   
306. bool LineIntersect(Line L1, Line L2) // 判断二维的两直线是否相交   
307. {   
308.     return (! ZERO((L1.p1 - L1.p2)*(L2.p1 - L2.p2)) );  // 是否平行   
309. }   
310. bool LineIntersect(Line3D L1, Line3D L2) // 判断三维的两直线是否相交   
311. {   
312.     Point3D p1 = L1.p1 - L1.p2;   
313.     Point3D p2 = L2.p1 - L2.p2;   
314.     Point3D p  = p1 * p2;   
315.     if (ZERO(p)) return false;      // 是否平行   
316.     p = (L2.p1 - L1.p2) * (L1.p1 - L1.p2);   
317.     return ZERO(p & L2.p2);         // 是否共面   
318. }   
319. bool LineSegIntersect(Line L1, Line L2) // 判断二维的两条线段是否相交   
320. {   
321.     return ( GEQ( max(L1.p1.x, L1.p2.x), min(L2.p1.x, L2.p2.x) ) &&   
322.         GEQ( max(L2.p1.x, L2.p2.x), min(L1.p1.x, L1.p2.x) ) &&   
323.         GEQ( max(L1.p1.y, L1.p2.y), min(L2.p1.y, L2.p2.y) ) &&   
324.         GEQ( max(L2.p1.y, L2.p2.y), min(L1.p1.y, L1.p2.y) ) &&   
325.         LEQ( ((L2.p1 - L1.p1) * (L1.p2 - L1.p1)) * ((L2.p2 -  L1.p1) * (L1.p2 - L1.p1)), 0 ) &&   
326.         LEQ( ((L1.p1 - L2.p1) * (L2.p2 - L2.p1)) * ((L1.p2 -  L2.p1) * (L2.p2 - L2.p1)), 0 ) );                
327. }   
328. bool LineSegIntersect(Line3D L1, Line3D L2) // 判断三维的两条线段是否相交   
329. {   
330.     // todo   
331.     return true;   
332. }   
333. // 计算两条二维直线的交点，结果在参数P中返回   
334. // 返回值说明了两条直线的位置关系:  COLINE   -- 共线  PARALLEL -- 平行  CROSS    -- 相交   
335. int CalCrossPoint(Line L1, Line L2, Point& P)   
336. {   
337.     double A1, B1, C1, A2, B2, C2;   
338.   
339.   
340.     A1 = L1.p2.y - L1.p1.y;   
341.     B1 = L1.p1.x - L1.p2.x;   
342.     C1 = L1.p2.x * L1.p1.y - L1.p1.x * L1.p2.y;   
343.   
344.   
345.     A2 = L2.p2.y - L2.p1.y;   
346.     B2 = L2.p1.x - L2.p2.x;   
347.     C2 = L2.p2.x * L2.p1.y - L2.p1.x * L2.p2.y;   
348.   
349.   
350.     if (EQ(A1 * B2, B1 * A2))    {   
351.         if (EQ( (A1 + B1) * C2, (A2 + B2) * C1 )) {   
352.             return COLINE;   
353.         } else {   
354.             return PARALLEL;   
355.         }   
356.     } else {   
357.         P.x = (B2 * C1 - B1 * C2) / (A2 * B1 - A1 * B2);   
358.         P.y = (A1 * C2 - A2 * C1) / (A2 * B1 - A1 * B2);   
359.         return CROSS;   
360.     }   
361. }   
362. // 计算两条三维直线的交点，结果在参数P中返回   
363. // 返回值说明了两条直线的位置关系 COLINE   -- 共线  PARALLEL -- 平行  CROSS    -- 相交  NONCOPLANAR -- 不公面   
364. int CalCrossPoint(Line3D L1, Line3D L2, Point3D& P)   
365. {   
366.     // todo   
367.     return 0;   
368. }   
369. // 计算点P到直线L的最近点   
370. Point NearestPointToLine(Point P, Line L)    
371. {   
372.     Point result;   
373.     double a, b, t;   
374.   
375.   
376.     a = L.p2.x - L.p1.x;   
377.     b = L.p2.y - L.p1.y;   
378.     t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b);   
379.   
380.   
381.     result.x = L.p1.x + a * t;   
382.     result.y = L.p1.y + b * t;   
383.     return result;   
384. }   
385. // 计算点P到线段L的最近点   
386. Point NearestPointToLineSeg(Point P, Line L)    
387. {   
388.     Point result;   
389.     double a, b, t;   
390.   
391.   
392.     a = L.p2.x - L.p1.x;   
393.     b = L.p2.y - L.p1.y;   
394.     t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b);   
395.   
396.   
397.     if ( GEQ(t, 0) && LEQ(t, 1) ) {   
398.         result.x = L.p1.x + a * t;   
399.         result.y = L.p1.y + b * t;   
400.     } else {   
401.         if ( Norm(P - L.p1) < Norm(P - L.p2) ) {   
402.             result = L.p1;   
403.         } else {   
404.             result = L.p2;   
405.         }   
406.     }   
407.     return result;   
408. }   
409. // 计算险段L1到线段L2的最短距离   
410. double MinDistance(Line L1, Line L2)    
411. {   
412.     double d1, d2, d3, d4;   
413.   
414.   
415.     if (LineSegIntersect(L1, L2)) {   
416.         return 0;   
417.     } else {   
418.         d1 = Norm( NearestPointToLineSeg(L1.p1, L2) - L1.p1 );   
419.         d2 = Norm( NearestPointToLineSeg(L1.p2, L2) - L1.p2 );   
420.         d3 = Norm( NearestPointToLineSeg(L2.p1, L1) - L2.p1 );   
421.         d4 = Norm( NearestPointToLineSeg(L2.p2, L1) - L2.p2 );   
422.            
423.         return min( min(d1, d2), min(d3, d4) );   
424.     }   
425. }   
426. // 求二维两直线的夹角，   
427. // 返回值是0~Pi之间的弧度   
428. double Inclination(Line L1, Line L2)   
429. {   
430.     Point u = L1.p2 - L1.p1;   
431.     Point v = L2.p2 - L2.p1;   
432.     return acos( (u & v) / (Norm(u)*Norm(v)) );   
433. }   
434. // 求三维两直线的夹角，   
435. // 返回值是0~Pi之间的弧度   
436. double Inclination(Line3D L1, Line3D L2)   
437. {   
438.     Point3D u = L1.p2 - L1.p1;   
439.     Point3D v = L2.p2 - L2.p1;   
440.     return acos( (u & v) / (Norm(u)*Norm(v)) );   
441. }   
442. /  
443.   
444.   
445. /  
446. // 判断两个矩形是否相交   
447. // 如果相邻不算相交   
448. bool Intersect(Rect_2 r1, Rect_2 r2)   
449. {   
450.     return ( max(r1.xl, r2.xl) < min(r1.xh, r2.xh) &&   
451.              max(r1.yl, r2.yl) < min(r1.yh, r2.yh) );   
452. }   
453. // 判断矩形r2是否可以放置在矩形r1内   
454. // r2可以任意地旋转   
455. //发现原来的给出的方法过不了OJ上的无归之室这题，  
456. //所以用了自己的代码  
457. bool IsContain(Rect r1, Rect r2)      //矩形的w>h  
458.  {   
459.      if(r1.w >r2.w && r1.h > r2.h) return true;  
460.      else  
461.      {  
462.         double r = sqrt(r2.w*r2.w + r2.h*r2.h) / 2.0;  
463.         double alpha = atan2(r2.h,r2.w);  
464.         double sita = asin((r1.h/2.0)/r);  
465.         double x = r * cos(sita - 2*alpha);  
466.         double y = r * sin(sita - 2*alpha);  
467.         if(x < r1.w/2.0 && y < r1.h/2.0 && x > 0 && y > -r1.h/2.0) return true;  
468.         else return false;  
469.      }  
470. }   
471.   
472.   
473.   
474.   
475. //圆  
476. Point Center(const Circle & C) //圆心  
477. {        
478.     return C.c;        
479. }      
480.   
481.   
482. double Area(const Circle &C)  
483. {  
484.  return pi*C.r*C.r;   
485. }   
486.   
487.   
488. double CommonArea(const Circle & A, const Circle & B) //两个圆的公共面积         
489. {        
490.     double area = 0.0;        
491.     const Circle & M = (A.r > B.r) ? A : B;        
492.     const Circle & N = (A.r > B.r) ? B : A;        
493.     double D = Distance(Center(M), Center(N));        
494.     if ((D < M.r + N.r) && (D > M.r - N.r))        
495.     {        
496.         double cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);        
497.         double cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);        
498.         double alpha = 2.0 * acos(cosM);        
499.         double beta  = 2.0 * acos(cosN);        
500.         double TM = 0.5 * M.r * M.r * sin(alpha);        
501.         double TN = 0.5 * N.r * N.r * sin(beta);        
502.         double FM = (alpha / (2*pi)) * Area(M);        
503.         double FN = (beta / (2*pi)) * Area(N);        
504.         area = FM + FN - TM - TN;        
505.     }        
506.     else if (D <= M.r - N.r)        
507.     {        
508.         area = Area(N);        
509.     }        
510.     return area;        
511. }   
512.        
513. bool IsInCircle(const Circle & C, const Rect_2 & rect)//判断圆是否在矩形内(不允许相切)  
514. {        
515.     return (GT(C.c.x - C.r, rect.xl)  
516.   &&  LT(C.c.x + C.r, rect.xh)  
517.   &&  GT(C.c.y - C.r, rect.yl)  
518.   &&  LT(C.c.y + C.r, rect.yh));        
519. }    
520.   
521.   
522. //判断2圆的位置关系  
523. //返回:   
524. //BAOHAN   = 1;        // 大圆包含小圆  
525. //NEIQIE   = 2;        // 内切  
526. //XIANJIAO = 3;        // 相交  
527. //WAIQIE   = 4;        // 外切  
528. //XIANLI   = 5;        // 相离  
529. int CirCir(const Circle &c1, const Circle &c2)//判断2圆的位置关系  
530. {  
531.  double dis = Distance(c1.c,c2.c);  
532.  if(LT(dis,fabs(c1.r-c2.r))) return BAOHAN;  
533.  if(EQ(dis,fabs(c1.r-c2.r))) return NEIQIE;  
534.  if(LT(dis,c1.r+c2.r) && GT(dis,fabs(c1.r-c2.r))) return XIANJIAO;  
535.  if(EQ(dis,c1.r+c2.r)) return WAIQIE;  
536.  return XIANLI;  
537. }  
538.   
539.   
540.   
541. int main()  
542. {  
543.  return 0;  
544. }  

1.14结构体表示几何图形

[cpp]  view plain  copy

1. //计算几何(二维)     
2. #include <cmath>     
3. #include <cstdio>     
4. #include <algorithm>     
5. using namespace std;     
6.   
7.   
8. typedef double TYPE;     
9. #define Abs(x) (((x)>0)?(x):(-(x)))     
10. #define Sgn(x) (((x)<0)?(-1):(1))     
11. #define Max(a,b) (((a)>(b))?(a):(b))     
12. #define Min(a,b) (((a)<(b))?(a):(b))     
13. #define Epsilon 1e-8     
14. #define Infinity 1e+10     
15. #define PI acos(-1.0)//3.14159265358979323846     
16. TYPE Deg2Rad(TYPE deg){return (deg * PI / 180.0);}     
17. TYPE Rad2Deg(TYPE rad){return (rad * 180.0 / PI);}     
18. TYPE Sin(TYPE deg){return sin(Deg2Rad(deg));}     
19. TYPE Cos(TYPE deg){return cos(Deg2Rad(deg));}     
20. TYPE ArcSin(TYPE val){return Rad2Deg(asin(val));}     
21. TYPE ArcCos(TYPE val){return Rad2Deg(acos(val));}     
22. TYPE Sqrt(TYPE val){return sqrt(val);}    
23.   
24.   
25. //点     
26. struct POINT     
27. {     
28.   TYPE x;     
29.   TYPE y;     
30.   POINT() : x(0), y(0) {};     
31.   POINT(TYPE _x_, TYPE _y_) : x(_x_), y(_y_) {};     
32. };     
33. // 两个点的距离     
34. TYPE Distance(const POINT & a, const POINT & b)     
35. {     
36.   return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));     
37. }     
38. //线段     
39. struct SEG     
40. {       
41.   POINT a; //起点     
42.   POINT b; //终点     
43.   SEG() {};     
44.   SEG(POINT _a_, POINT _b_):a(_a_),b(_b_) {};     
45. };       
46. //直线(两点式)     
47. struct LINE     
48. {     
49.   POINT a;     
50.   POINT b;     
51.   LINE() {};     
52.   LINE(POINT _a_, POINT _b_) : a(_a_), b(_b_) {};     
53. };     
54. //直线(一般式)     
55. struct LINE2     
56. {     
57.   TYPE A,B,C;     
58.   LINE2() {};     
59.   LINE2(TYPE _A_, TYPE _B_, TYPE _C_) : A(_A_), B(_B_), C(_C_) {};     
60. };     
61.   
62.   
63. //两点式化一般式     
64. LINE2 Line2line(const LINE & L) // y=kx+c k=y/x  
65. {     
66.   LINE2 L2;     
67.   L2.A = L.b.y - L.a.y;     
68.   L2.B = L.a.x - L.b.x;     
69.   L2.C = L.b.x * L.a.y - L.a.x * L.b.y;     
70.   return L2;     
71. }     
72.   
73.   
74. // 引用返回直线 Ax + By + C =0 的系数     
75. void Coefficient(const LINE & L, TYPE & A, TYPE & B, TYPE & C)     
76. {     
77.   A = L.b.y - L.a.y;     
78.   B = L.a.x - L.b.x;     
79.   C = L.b.x * L.a.y - L.a.x * L.b.y;     
80. }     
81. void Coefficient(const POINT & p,const TYPE a,TYPE & A,TYPE & B,TYPE & C)     
82. {     
83.   A = Cos(a);     
84.   B = Sin(a);     
85.   C = - (p.y * B + p.x * A);     
86. }     
87. /判等(值，点，直线)     
88. bool IsEqual(TYPE a, TYPE b)     
89. {     
90.   return (Abs(a - b) <Epsilon);     
91. }     
92. bool IsEqual(const POINT & a, const POINT & b)     
93. {     
94.   return (IsEqual(a.x, b.x) && IsEqual(a.y, b.y));     
95. }     
96. bool IsEqual(const LINE & A, const LINE & B)     
97. {     
98.   TYPE A1, B1, C1;     
99.   TYPE A2, B2, C2;     
100.   Coefficient(A, A1, B1, C1);     
101.   Coefficient(B, A2, B2, C2);     
102.   return IsEqual(A1 * B2, A2 * B1) && IsEqual(A1 * C2, A2 * C1) && IsEqual(B1 * C2, B2 * C1);     
103. }     
104. // 矩形     
105. struct RECT     
106. {     
107.   POINT a; // 左下点       
108.   POINT b; // 右上点       
109.   RECT() {};     
110.   RECT(const POINT & _a_, const POINT & _b_) { a = _a_; b = _b_; }     
111. };     
112.   
113.   
114. //矩形化标准     
115. RECT Stdrect(const RECT & q)  
116. {     
117.   TYPE t;     
118.   RECT p=q;     
119.   if(p.a.x > p.b.x) swap(p.a.x , p.b.x);      
120.   if(p.a.y > p.b.y) swap(p.a.y , p.b.y);      
121.   return p;     
122. }     
123.   
124.   
125. //根据下标返回矩形的边       
126. SEG Edge(const RECT & rect, int idx)     
127. {     
128.   SEG edge;     
129.   while (idx < 0) idx += 4;     
130.   switch (idx % 4)     
131.   {     
132.   case 0: //下边  
133.     edge.a = rect.a;     
134.     edge.b = POINT(rect.b.x, rect.a.y);     
135.     break;     
136.   case 1: //右边  
137.     edge.a = POINT(rect.b.x, rect.a.y);     
138.     edge.b = rect.b;     
139.     break;     
140.   case 2: //上边    
141.     edge.a = rect.b;     
142.     edge.b = POINT(rect.a.x, rect.b.y);     
143.     break;     
144.   case 3: //左边    
145.     edge.a = POINT(rect.a.x, rect.b.y);     
146.     edge.b = rect.a;     
147.     break;     
148.   default:     
149.     break;     
150.   }     
151.   return edge;     
152. }     
153.   
154.   
155. //矩形的面积     
156. TYPE Area(const RECT & rect)     
157. {     
158.   return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y);     
159. }     
160.   
161.   
162. //两个矩形的公共面积       
163. TYPE CommonArea(const RECT & A, const RECT & B)     
164. {     
165.   TYPE area = 0.0;     
166.   POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y));     
167.   POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y));     
168.   if( (LL.x <= UR.x) && (LL.y <= UR.y) )     
169.   {     
170.     area = Area(RECT(LL, UR));     
171.   }     
172.   return area;     
173. }    
174. //判断圆是否在矩形内(不允许相切)     
175. bool IsInCircle(const CIRCLE & circle, const RECT & rect)     
176. {     
177.   return (circle.x - circle.r > rect.a.x) &&     
178.     (circle.x + circle.r < rect.b.x) &&     
179.     (circle.y - circle.r > rect.a.y) &&     
180.     (circle.y + circle.r < rect.b.y);     
181. }     
182.   
183.   
184. //判断矩形是否在圆内(不允许相切)     
185. bool IsInRect(const CIRCLE & circle, const RECT & rect)     
186. {     
187.   POINT c,d;     
188.   c.x=rect.a.x; c.y=rect.b.y;     
189.   d.x=rect.b.x; d.y=rect.a.y;     
190.   return (Distance( Center(circle) , rect.a ) < circle.r) &&     
191.     (Distance( Center(circle) , rect.b ) < circle.r) &&     
192.     (Distance( Center(circle) , c ) < circle.r) &&     
193.     (Distance( Center(circle) , d ) < circle.r);     
194. }     
195.   
196.   
197. //判断矩形是否与圆相离(不允许相切)     
198. bool Isoutside(const CIRCLE & circle, const RECT & rect)     
199. {     
200.   POINT c,d;     
201.   c.x=rect.a.x; c.y=rect.b.y;     
202.   d.x=rect.b.x; d.y=rect.a.y;     
203.   return (Distance( Center(circle) , rect.a ) > circle.r) &&     
204.     (Distance( Center(circle) , rect.b ) > circle.r) &&     
205.     (Distance( Center(circle) , c ) > circle.r) &&     
206.     (Distance( Center(circle) , d ) > circle.r) &&     
207.     (rect.a.x > circle.x || circle.x > rect.b.x || rect.a.y > circle.y || circle.y > rect.b.y) ||     
208.     ((circle.x - circle.r > rect.b.x) ||     
209.     (circle.x + circle.r < rect.a.x) ||     
210.     (circle.y - circle.r > rect.b.y) ||     
211.     (circle.y + circle.r < rect.a.y));     
212. }     

1.15四城部分几何模板

[cpp]  view plain  copy

1. /* 
2. 1.注意实际运用的时候可以用sqrd代替dist提高精度，节省时间 
3. */  
4. #include <iostream>  
5. #include <math.h>  
6. #include <algorithm>  
7. using namespace std;  
8.   
9.   
10. const double INF = 10e300;  
11. const double EPS = 1e-8;  
12. const double PI = acos(-1.0);  
13.   
14.   
15. inline int dblcmp(double a, double b) {if(fabs(a-b) < EPS) return 0;if(a < b) return -1;return 1;}  
16. inline double Max(double a, double b) { if(dblcmp(a, b) == 1) return a; return b; }  
17. inline double Min(double a, double b) { if(dblcmp(a, b) == 1) return b; return a; }  
18. inline double Agl(double deg) { return deg * PI / 180.0; }  
19.   
20.   
21. struct Point { double x, y; void set(double a, double b) { x = a; y = b; } };  
22. struct Vec { double x, y; void set(Point& a, Point& b) { x = b.x-a.x; y = b.y-a.y; } };  
23. struct Line { double a, b, c; Point st, end;  
24. void set(Point& u, Point& v) {a = v.y - u.y; b = u.x - v.x; c = a*u.x + b*u.y; st = u; end = v; } };  
25.   
26.   
27. inline double dist(Point& a, Point& b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); }  
28. inline double sqrd(Point& a, Point& b) { return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y); }  
29. inline double dot(Vec& a, Vec& b) { return a.x * b.x + a.y * b.y; }  
30. inline double cross(Vec& a, Vec& b) { return a.x * b.y - a.y * b.x; }  
31. inline double cross(Point& a, Point& b, Point& c) {Vec x, y; x.set(a, b); y.set(a, c); return cross(x, y); }  
32. //返回1代表a在bc之间 0代表在端点 -1代表在外面  
33. inline int between(Point& a, Point& b, Point& c) { Vec x, y; x.set(a,b); y.set(a,c); return dblcmp(dot(x, y),0); }  
34.   
35.   
36. //3维坐标转换 输入是度数  
37. void trans(double lat, double log, double& x, double& y, double& z, double radius) {  
38.     x = radius * cos(lat) * cos(log);  
39.     y = radius * cos(lat) * sin(log);  
40.     z = radius * sin(lat);  
41. }  
42.   
43.   
44. //求两点的平分线  
45. Line bisector(Point& a, Point& b) {  
46.     Line ab, ans; ab.set(a, b);  
47.     double midx = (a.x + b.x)/2.0,  midy = (a.y + b.y)/2.0;  
48.     ans.a = -ab.b, ans.b = -ab.a, ans.c = -ab.b * midx + ab.a * midy;  
49.     return ans;  
50. }  
51.   
52.   
53.   
54.   
55. //线线相交 如果平行 返回-1, 重合返回-2  
56. int line_line_intersect(Line& l1, Line& l2, Point& s) {  
57.     double det = l1.a*l2.b - l2.a*l1.b;  
58.     if(dblcmp(det, 0.0) == 0) { //平行或者重合  
59.         if(dblcmp(point_line_dist(l1.st, l2.st, l2.end, 0), 0) == 0)   
60.             return -2;  
61.         return -1;  
62.     }  
63.     s.x = (l2.b*l1.c - l1.b*l2.c)/det;  
64.     s.y = (l1.a*l2.c - l2.a*l1.c)/det;  
65.     return 1;  
66. }  
67.   
68.   
69. //2线段相交 ab, cd 交点是s 平行返回-1, 重合返回-2, 不在线段上面返回0 在线段中间返回1 在线段两端返回2  
70. int seg_seg_intersect(Point& a, Point& b, Point& c, Point& d, Point& s) {  
71.     Line l1, l2; l1.set(a, b); l2.set(c, d);  
72.     int ans = line_line_intersect(l1, l2, s);  
73.     if(ans == 1) {  
74.         if(between(s, a, b) == 1 && between(s, c, d) == 1)   
75.             return 1;  
76.         if(between(s, a, b) == -1 && between(s, c, d) == -1)  
77.             return 0;  
78.         return 2;  
79.     }  
80.     return ans;  
81. }  
82.   
83.   
84. //求三点共圆 中心放在center中 返回半径  
85. double center_3point(Point& a, Point& b, Point& c, Point& center) {  
86.     Line x = bisector(a, b), y = bisector(b, c);  
87.     line_line_intersect(x, y, center);  
88.     return dist(center, a);  
89. }  

1.16  一些代码

1.16.1 最小圆覆盖_zju1450

/*
包含点集所有点的最小圆的算法
最小圆覆盖 


 http://acm.zju.edu.cn/show_problem.php?pid=1450
 
 相关题目最小球包含  http://acm.pku.edu.cn/JudgeOnline/problem?id=2069 


平面上有n个点，给定n个点的坐标，试找一个半径最小的圆，将n
个点全部包围，点可以在圆上。 


1. 在点集中任取3点A,B,C。
2. 作一个包含A,B,C三点的最小圆,圆周可能通过这3点，也可能只通过
其中两点,但包含第3点.后一种情况圆周上的两点一定是位于圆的一条直
径的两端。
3. 在点集中找出距离第2步所建圆圆心最远的D点，若D点已在圆内或圆周上，
则该圆即为所求的圆，算法结束.则，执行第4步。
4. 在A,B,C,D中选3个点,使由它们生成的一个包含这4个点的圆为最小，这3
点成为新的A,B,C，返回执行第2步。若在第4步生成的圆的圆周只通过A,B,C,D
中的两点，则圆周上的两点取成新的A和B,从另两点中任取一点作为新的C。


程序设计题解上的解题报告：
对于一个给定的点集A，记MinCircle(A)为点集A的最小外接圆，显然，对于所
有的点集情况A,MinCircle(A)都是存在且惟一的。需要特别说明的是，当A为空
集时，MinCircle(A)为空集，当A={a}时，MinCircle(A)圆心坐标为a，半径为0；
   显然，MinCircle(A)可以有A边界上最多三个点确定(当点集A中点的个数大于
1时，有可能两个点确定了MinCircle(A))，也就是说存在着一个点集B，|B|<=3
且B包含与A，有MinCircle(B)=MinCircle(A).所以，如果a不属于B，则
MinCircle(A-{a})=MinCircle(A);如果MinCircle(A-{a})不等于MinCircle(A),则
a属于B。
    所以我们可以从一个空集R开始，不断的把题目中给定的点集中的点加入R，同
时维护R的外接圆最小，这样就可以得到解决该题的算法。 

pku2069

[cpp]  view plain  copy

1. */  
2. #include <stdio.h>  
3. #include <math.h>  
4.   
5.   
6. const int maxn = 1005;  
7. //const double eps = 1e-6;  
8.   
9.   
10. struct TPoint  
11. {  
12.     double x, y;  
13.     TPoint operator-(TPoint &a)  
14.     {  
15.         TPoint p1;  
16.         p1.x = x - a.x;  
17.         p1.y = y - a.y;  
18.         return p1;  
19.     }  
20. };  
21.   
22.   
23. struct TCircle  
24. {  
25.     double r;  
26.     TPoint centre;  
27. };  
28.   
29.   
30. struct TTriangle  
31. {  
32.     TPoint t[3];  
33. };  
34.   
35.   
36. TCircle c;  
37. TPoint a[maxn];  
38.   
39.   
40. double distance(TPoint p1, TPoint p2)  
41. {  
42.     TPoint p3;  
43.     p3.x = p2.x - p1.x;  
44.     p3.y = p2.y - p1.y;  
45.     return sqrt(p3.x * p3.x + p3.y * p3.y);  
46. }  
47.   
48.   
49. double triangleArea(TTriangle t)  
50. {  
51.     TPoint p1, p2;  
52.     p1 = t.t[1] - t.t[0];  
53.     p2 = t.t[2] - t.t[0];  
54.     return fabs(p1.x * p2.y - p1.y * p2.x) / 2;  
55. }  
56.   
57.   
58. TCircle circumcircleOfTriangle(TTriangle t)  
59. {  
60.     //三角形的外接圆  
61.     TCircle tmp;  
62.     double a, b, c, c1, c2;  
63.     double xA, yA, xB, yB, xC, yC;  
64.     a = distance(t.t[0], t.t[1]);  
65.     b = distance(t.t[1], t.t[2]);  
66.     c = distance(t.t[2], t.t[0]);  
67.     //根据S = a * b * c / R / 4;求半径R   
68.     tmp.r = a * b * c / triangleArea(t) / 4;  
69.       
70.     xA = t.t[0].x; yA = t.t[0].y;  
71.     xB = t.t[1].x; yB = t.t[1].y;  
72.     xC = t.t[2].x; yC = t.t[2].y;  
73.     c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;  
74.     c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;  
75.       
76.     tmp.centre.x = (c1 * (yA - yC) - c2 * (yA - yB)) /   
77.          ((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));   
78.     tmp.centre.y = (c1 * (xA - xC) - c2 * (xA - xB)) /   
79.          ((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));   
80.            
81.     return tmp;       
82. }  
83.   
84.   
85. TCircle MinCircle2(int tce, TTriangle ce)  
86. {  
87.     TCircle tmp;  
88.     if(tce == 0) tmp.r = -2;  
89.     else if(tce == 1)   
90.     {  
91.         tmp.centre = ce.t[0];  
92.         tmp.r = 0;  
93.     }  
94.     else if(tce == 2)  
95.     {  
96.         tmp.r = distance(ce.t[0], ce.t[1]) / 2;  
97.         tmp.centre.x = (ce.t[0].x + ce.t[1].x) / 2;  
98.         tmp.centre.y = (ce.t[0].y + ce.t[1].y) / 2;       
99.     }  
100.     else if(tce == 3) tmp = circumcircleOfTriangle(ce);  
101.     return tmp;  
102. }  
103.   
104.   
105. void MinCircle(int t, int tce, TTriangle ce)  
106. {  
107.     int i, j;  
108.     TPoint tmp;  
109.     c = MinCircle2(tce, ce);  
110.     if(tce == 3) return;  
111.     for(i = 1;i <= t;i++)  
112.     {  
113.         if(distance(a[i], c.centre) > c.r)  
114.         {  
115.             ce.t[tce] = a[i];  
116.             MinCircle(i - 1, tce + 1, ce);  
117.             tmp = a[i];  
118.             for(j = i;j >= 2;j--)  
119.             {  
120.                 a[j] = a[j - 1];  
121.             }  
122.             a[1] = tmp;  
123.         }  
124.     }  
125. }  
126.   
127.   
128. void run(int n)  
129. {  
130.     TTriangle ce;  
131.     int i;  
132.     MinCircle(n, 0, ce);  
133.     printf("%.2lf %.2lf %.2lf\n", c.centre.x, c.centre.y, c.r);  
134. }  
135.   
136.   
137. int main()  
138. {  
139.     freopen("circle.in", "r", stdin);  
140.     freopen("out.txt", "w", stdout);  
141.     int n;  
142.     while(scanf("%d", &n) != EOF && n)  
143.     {  
144.         for(int i = 1;i <= n;i++)  
145.             scanf("%lf%lf", &a[i].x, &a[i].y);  
146.         run(n);  
147.     }  
148.     return 0;  
149. }  

1.16.2 直线旋转_两凸包的最短距离(poj3608)

[cpp]  view plain  copy

1. #include <stdio.h>  
2. #include <math.h>  
3.   
4.   
5. #define pi acos(-1.0)  
6. #define eps 1e-6  
7. #define inf 1e250  
8. #define Maxn 10005  
9.   
10.   
11. typedef struct TPoint  
12. {  
13.     double x, y;  
14. }TPoint;  
15.   
16.   
17. typedef struct TPolygon  
18. {  
19.     TPoint p[Maxn];  
20.     int n;  
21. }TPolygon;  
22.   
23.   
24. typedef struct TLine  
25. {  
26.     double a, b, c;  
27. }TLine;  
28.   
29.   
30. double max(double a, double b)  
31. {  
32.     if(a > b) return a;  
33.     return b;  
34. }  
35.   
36.   
37. double min(double a, double b)  
38. {  
39.     if(a < b) return a;  
40.     return b;  
41. }  
42.   
43.   
44. double distance(TPoint p1, TPoint p2)  
45. {  
46.     return sqrt((p1.x - p2.x) * (p1.x - p2.x)   
47.      + (p1.y - p2.y) * (p1.y - p2.y));  
48. }  
49.   
50.   
51. TLine lineFromSegment(TPoint p1, TPoint p2)  
52. {  
53.     TLine tmp;  
54.     tmp.a = p2.y - p1.y;  
55.     tmp.b = p1.x - p2.x;  
56.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
57.     return tmp;  
58. }  
59.   
60.   
61. double polygonArea(TPolygon p)  
62. {  
63.     int i, n;  
64.     double area;  
65.     n = p.n;  
66.     area = 0;  
67.     for(i = 1;i <= n;i++)  
68.         area += (p.p[i - 1].x * p.p[i % n].y - p.p[i % n].x * p.p[i - 1].y);  
69.   
70.   
71.     return area / 2;    
72. }  
73.   
74.   
75. void ChangeClockwise(TPolygon &polygon)  
76. {  
77.     TPoint tmp;  
78.     int i;  
79.     for(i = 0;i <= (polygon.n - 1) / 2;i++)  
80.     {  
81.         tmp = polygon.p[i];  
82.         polygon.p[i] = polygon.p[polygon.n - 1 - i];  
83.         polygon.p[polygon.n - 1 - i] = tmp;           
84.     }  
85. }  
86.   
87.   
88. double disPointToSeg(TPoint p1, TPoint p2, TPoint p3)  
89. {  
90.     double a = distance(p1, p2);  
91.     double b = distance(p1, p3);  
92.     double c = distance(p2, p3);  
93.     if(fabs(a + b - c) < eps) return 0;  
94.     if(fabs(a + c - b) < eps || fabs(b + c - a) < eps) return min(a, b);  
95.     double t1 = -a * a + b * b + c * c;  
96.     double t2 = a * a - b * b + c * c;  
97.     if(t1 <= 0 || t2 <= 0) return min(a, b);  
98.       
99.     TLine l1 = lineFromSegment(p2, p3);  
100.     return fabs(l1.a * p1.x + l1.b * p1.y + l1.c) / sqrt(l1.a * l1.a + l1.b * l1.b);        
101. }  
102.   
103.   
104. double disPallSeg(TPoint p1, TPoint p2, TPoint p3, TPoint p4)  
105. {  
106.     return min(min(disPointToSeg(p1, p3, p4), disPointToSeg(p2, p3, p4)),  
107.      min(disPointToSeg(p3, p1, p2), disPointToSeg(p4, p1, p2)));  
108. }  
109.   
110.   
111. double angle(TPoint p1, TPoint p2, double SlewRate)  
112. {  
113.     double ang, tmp;  
114.     TPoint p;  
115.     p.x = p2.x - p1.x;  
116.     p.y = p2.y - p1.y;  
117.     if(fabs(p.x) < eps)  
118.     {  
119.         if(p.y > 0) ang = pi / 2;  
120.         else ang = 3 * pi / 2;  
121.     }  
122.     else   
123.     {  
124.         ang = atan(p.y / p.x);  
125.         if(p.x < 0) ang += pi;  
126.     }  
127.     while(ang < 0) ang += 2 * pi;  
128.     if(ang >= pi) SlewRate += pi;  
129.     if(ang > SlewRate) tmp = ang - SlewRate;  
130.     else tmp = pi - (SlewRate - ang);  
131.     while(tmp >= pi) tmp -= pi;  
132.     if(fabs(tmp - pi) < eps) tmp = 0;  
133.     return tmp;  
134. }  
135.   
136.   
137. int main()  
138. {  
139.     int n, m, i;  
140.     TPolygon polygon1, polygon2;  
141.     double ymin1, ymax2, ans, d;  
142.     int k1, k2;  
143.     while(scanf("%d%d", &n, &m) && n)  
144.     {  
145.         polygon1.n = n;  
146.         polygon2.n = m;  
147.         for(i = 0;i < n;i++)  
148.             scanf("%lf%lf", &polygon1.p[i].x, &polygon1.p[i].y);  
149.         for(i = 0;i < m;i++)  
150.             scanf("%lf%lf", &polygon2.p[i].x, &polygon2.p[i].y);      
151.         if(polygonArea(polygon1) < 0) ChangeClockwise(polygon1);  
152.         if(polygonArea(polygon2) < 0) ChangeClockwise(polygon2);  
153.         ymin1 = inf, ymax2 = -inf;  
154.         for(i = 0;i < n;i++)  
155.             if(polygon1.p[i].y < ymin1) ymin1 = polygon1.p[i].y , k1 = i;  
156.         for(i = 0;i < m;i++)  
157.             if(polygon2.p[i].y > ymax2) ymax2 = polygon2.p[i].y , k2 = i;      
158.         double SlewRate = 0;  
159.         double angle1, angle2;  
160.         ans = inf;  
161.         double Slope = 0;  
162.         while(Slope <= 360)    
163.         {     
164.             while(SlewRate >= pi) SlewRate -= pi;  
165.             if(fabs(pi - SlewRate) < eps) SlewRate = 0;  
166.             angle1 = angle(polygon1.p[k1], polygon1.p[(k1 + 1) % n], SlewRate);  
167.             angle2 = angle(polygon2.p[k2], polygon2.p[(k2 + 1) % m], SlewRate);   
168.             if(fabs(angle1 - angle2) < eps)  
169.             {  
170.                 d = disPallSeg(polygon1.p[k1], polygon1.p[(k1 + 1) % n], polygon2.p[k2], polygon2.p[(k2 + 1) % m]);   
171.                 if(d < ans) ans = d;  
172.                 k1++;  
173.                 k1 %= n;  
174.                 k2++;  
175.                 k2 %= m;   
176.                 SlewRate += angle1;  
177.                 Slope += angle1;  
178.             }  
179.             else if(angle1 < angle2)  
180.             {  
181.                 d = disPointToSeg(polygon2.p[k2], polygon1.p[k1], polygon1.p[(k1 + 1) % n]);  
182.                 if(d < ans) ans = d;  
183.                 k1++;  
184.                 k1 %= n;  
185.                 SlewRate += angle1;  
186.                 Slope += angle1;  
187.             }  
188.             else   
189.             {  
190.                 d = disPointToSeg(polygon1.p[k1], polygon2.p[k2], polygon2.p[(k2 + 1) % m]);  
191.                 if(d < ans) ans = d;  
192.                 k2++;  
193.                 k2 %= m;  
194.                 SlewRate += angle2;  
195.                 Slope += angle2;  
196.             }  
197.         }  
198.         printf("%.5lf\n", ans);  
199.     }  
200.     return 0;  
201. }  

1.16.3 扇形的重心

[cpp]  view plain  copy

1. //Xc = 2*R*sinA/3/A   
2. //A为圆心角的一半  
3. #include <stdio.h>  
4. #include <math.h>  
5. int main()  
6. {  
7.     double r, angle;  
8.     while(scanf("%lf%lf", &r, &angle) != EOF){  
9.         angle /= 2;  
10.         printf("%.6lf\n", 2 * r * sin(angle) / 3 / angle);  
11.     }  
12.     return 0;  
13. }  

1.16.4 根据经度纬度求球面距离

[cpp]  view plain  copy

1. /* 
2. 假设地球是球体，  
3. 设地球上某点的经度为lambda,纬度为phi，  
4. 则这点的空间坐标是  
5. x=cos(phi)*cos(lambda)  
6. y=cos(phi)*sin(lambda)  
7. z=sin(phi)  
8. 设地球上两点的空间坐标分别为(x1,y1,z1),(x2,y2,z2)  
9. 直线距离即为R*sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)), 
10. 则它们的夹角为  
11. A = acos(x1 * x2 + y1 * y2 + z1 * z2)， 
12. 则两地距离为 A * R，其中R为地球平均半径6371  
13. */  
14.   
15.   
16. /* 
17. 这里坐标都要乘以半径R，但由于是求角度，所以统一都没有乘  
18. 注意这里还要判断坐标的正负和经度纬度的规定有关  
19.  
20.  
21. pku_3407 
22. */  
23. #include <stdio.h>  
24. #include <math.h>  
25.   
26.   
27. const double pi = acos(-1.0);  
28.   
29.   
30. struct TPoint  
31. {  
32.    double x, y, z;  
33. };  
34.    
35. int main()  
36. {  
37.     double w1, wm1, j1, jm1, wd1, wd2;  
38.     double w2, wm2, j2, jm2, jd1, jd2;  
39.     TPoint p1, p2;  
40.     char chr1, chr2;  
41.     while(scanf("%lf%lf ", &w1, &wm1) != EOF){  
42.         scanf("%c ", &chr1);  
43.         scanf("%lf %lf %c", &j1, &jm1, &chr2);  
44.         wd1 = (w1 + wm1 / 60) * pi / 180;  
45.         jd1 = (j1 + jm1 / 60) * pi / 180;  
46.         if(chr1 == 'S') wd1 *= -1.0;  
47.         if(chr2 == 'W') jd1 *= -1.0;  
48.         p1.x = cos(wd1) * cos(jd1);  
49.         p1.y = cos(wd1) * sin(jd1);  
50.         p1.z = sin(wd1);  
51.         scanf("%lf %lf %c %lf %lf %c", &w2, &wm2, &chr1, &j2, &jm2, &chr2);  
52.         wd2 = (w2 + wm2 / 60) * pi / 180;  
53.         jd2 = (j2 + jm2 / 60) * pi / 180;  
54.         if(chr1 == 'S') wd2 *= -1.0;  
55.         if(chr2 == 'W') jd2 *= -1.0;  
56.         p2.x = cos(wd2) * cos(jd2);  
57.         p2.y = cos(wd2) * sin(jd2);  
58.         p2.z = sin(wd2);  
59.         double a = acos(p1.x * p2.x + p1.y * p2.y + p1.z * p2.z);  
60.         printf("%.3lf\n", a * 6370.0);     
61.     }  
62.     return 0;  
63. }  

1.16.5 多边形的重心

[cpp]  view plain  copy

1. /* 
2. 题目描述： 
3. 有一个密度均匀的平面N多边形(3 <= N <= 1000000)，可能凹也可能凸，但没有边相交叉， 
4. 另外已知N个有序(顺时针或逆时针)顶点的坐标值，第j个顶点坐标为（Xi ， Yi ），且满 
5. 足 (|Xi|, |Yi| <= 20000)，求这个平面多边形的重心。 
6.  
7.  
8.   解题过程： 
9.   从第1个顶点出发，分别连接第i, i+1个顶点组成三角形Ti,1 < i < n,  
10.   一共n-2个三角形正好是多连形的一个划分，分别求出每个三角形的面积Si， 
11.   总面积为各个面积相加 
12.   根据物理学知识得：n个点(xi,yi)每个重量是mi,则重心是    
13.   X = (x1*M1+x2*M2+...+xn*Mn) / (M1+M2+....+Mn)    
14.   Y = (y1*M1+y2*M2+...+yn*Mn) / (M1+M2+....+Mn)    
15.   另个需要用的知识有： 
16.   已知3点求三角形的面积，设三点分别为p[0].x, p[0].y p[1].x, p[1].y p[1].x, p[1].y 
17.   面积s =[ p[0].x*p[1].y-p[1].x*p[0].y + p[1].x*p[2].y-p[2].x*p[1].y + p[2].x*p[0].y-p[0].x*p[2].y ] / 2 ,  
18.   这是这3个点是逆时针的值，顺时针取负。 
19.   已知3点求重心，x = (p[0].x+p[1].x+p[2].x)/3.0 , y = (p[0].y+p[1].y+p[2].y)/3.0  
20.   另外在求解的过程中，不需要考虑点的输入顺序是顺时针还是逆时针，相除后就抵消了， 
21.   还要注意 一点是不必在求每个小三角形的重心时都除以3，可以在最后除一下 
22.      
23. */  
24. /*fzu_1132*/  
25. #include <stdio.h>  
26. #include <math.h>  
27.   
28.   
29. typedef struct TPoint  
30. {  
31.     double x;  
32.     double y;  
33. }TPoint;  
34.   
35.   
36. double triangleArea(TPoint p0, TPoint p1, TPoint p2)  
37. {  
38.     //已知三角形三个顶点的坐标，求三角形的面积   
39.     double k = p0.x * p1.y + p1.x * p2.y  
40.         + p2.x * p0.y - p1.x * p0.y   
41.         - p2.x * p1.y - p0.x * p2.y;  
42.     //if(k >= 0) return k / 2;  
43. //  else return -k / 2;      
44.        return k / 2;  
45. }  
46.   
47.   
48. int main()  
49. {  
50.     int i, n, test;  
51.     TPoint p0, p1, p2, center;  
52.     double area, sumarea, sumx, sumy;      
53.     scanf("%d", &test);    
54.     while(test--){  
55.         scanf("%d", &n);  
56.         scanf("%lf%lf", &p0.x, &p0.y);  
57.         scanf("%lf%lf", &p1.x, &p1.y);  
58.         sumx = 0;  
59.         sumy = 0;  
60.         sumarea = 0;  
61.         for(i = 2;i < n;i++){  
62.             scanf("%lf%lf", &p2.x, &p2.y);  
63.             center.x = p0.x + p1.x + p2.x;  
64.             center.y = p0.y + p1.y + p2.y;    
65.             area =  triangleArea(p0, p1, p2);  
66.             sumarea += area;  
67.             sumx += center.x * area;  
68.             sumy += center.y * area;   
69.             p1 = p2;             
70.         }  
71.         printf("%.2lf %.2lf\n", sumx / sumarea / 3, sumy / sumarea / 3);  
72.     }  
73.     return 0;  
74. }  

1.16.6 存不存在一个平面把两堆点分开(poj3643)

[cpp]  view plain  copy

1. #include <stdio.h>  
2. struct point  
3. {  
4.     double x, y, z;  
5. }pa[201], pb[201];  
6. int main()   
7. {   
8.     int n, m, i;   
9.     while (scanf("%d", &n), n != -1)   
10.     {   
11.         for (i = 0; i < n; i++)   
12.             scanf("%lf%lf%lf", &pa[i].x, &pa[i].y, &pa[i].z);   
13.         scanf("%d", &m);   
14.         for (i = 0; i < m; i++)   
15.             scanf("%lf%lf%lf", &pb[i].x, &pb[i].y, &pb[i].z);  
16.         int cnt = 0, finish = 0;   
17.         double a = 0, b = 0, c = 0, d = 0;   
18.         while (cnt < 100000 && !finish)  
19.         {   
20.             finish = 1;   
21.             for (i = 0; i < n; i++)   
22.                 if (a * pa[i].x + b * pa[i].y + c * pa[i].z + d > 0)   
23.                 {   
24.                     a -= pa[i].x;   
25.                     b -= pa[i].y;   
26.                     c -= pa[i].z;   
27.                     d -= 3;   
28.                     finish = 0;   
29.                 }  
30.             for (i = 0; i < m; i++)   
31.                 if (a * pb[i].x + b * pb[i].y + c * pb[i].z + d <= 0)   
32.                 {   
33.                     a += pb[i].x;   
34.                     b += pb[i].y;   
35.                     c += pb[i].z;   
36.                     d += 3;   
37.                     finish = 0;   
38.                 }  
39.             cnt++;   
40.         }  
41.         printf("%lf %lf %lf %lf\n", a, b, c, d);   
42.     }  
43.     return 0;  
44. }  

1.16.7 pku_3335_判断多边形的核是否存在

[cpp]  view plain  copy

1. /*多边形的核*/  
2. #include <stdio.h>  
3. #include <math.h>  
4.   
5.   
6. #define Maxn 3005  
7. const double eps = 1e-10;  
8.   
9.   
10. typedef struct TPodouble   
11. {  
12.     double x;  
13.     double y;  
14. }TPoint;  
15.   
16.   
17. typedef struct TPolygon  
18. {  
19.     TPoint p[Maxn];  
20.     int n;  
21. };  
22.   
23.   
24. typedef struct TLine  
25. {  
26.     double a, b, c;  
27. }TLine;  
28.   
29.   
30. bool same(TPoint p1, TPoint p2)  
31. {  
32.     if(p1.x != p2.x) return false;  
33.     if(p1.y != p2.y) return false;  
34.     return true;  
35. }  
36.    
37. double multi(TPoint p1, TPoint p2, TPoint p0)  
38. {  
39.     //求矢量[p0, p1], [p0, p2]的叉积  
40.     //p0是顶点   
41.     return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);  
42.     //若结果等于0，则这三点共线  
43.     //若结果大于0，则p0p2在p0p1的逆时针方向  
44.     //若结果小于0，则p0p2在p0p1的顺时针方向   
45. }  
46.   
47.   
48. TLine lineFromSegment(TPoint p1, TPoint p2)  
49. {  
50.     //线段所在直线,返回直线方程的三个系统   
51.     TLine tmp;  
52.     tmp.a = p2.y - p1.y;  
53.     tmp.b = p1.x - p2.x;  
54.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
55.     return tmp;  
56. }  
57.   
58.   
59. TPoint LineInter(TLine l1, TLine l2)  
60. {  
61.     //求两直线得交点坐标  
62.     TPoint tmp;   
63.     double a1 = l1.a;  
64.     double b1 = l1.b;  
65.     double c1 = l1.c;  
66.     double a2 = l2.a;  
67.     double b2 = l2.b;  
68.     double c2 = l2.c;  
69.     //注意这里b1 = 0   
70.     if(fabs(b1) < eps){  
71.         tmp.x = -c1 / a1;    
72.         tmp.y = (-c2 - a2 * tmp.x) / b2;  
73.     }         
74.     else{  
75.         tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);  
76.         tmp.y = (-c1 - a1 * tmp.x) / b1;  
77.     }  
78.     return tmp;  
79. }  
80.   
81.   
82. TPolygon Cut_polygon(TPoint p1, TPoint p2, TPolygon polygon)  
83. {  
84.     TPolygon new_polygon;  
85.     TPoint interp;  
86.     TLine l1, l2;  
87.     int i, j;  
88.     double t1, t2;  
89.     new_polygon.n = 0;  
90.     for(i = 0;i <= polygon.n - 1;i++){  
91.         t1 = multi(p2, polygon.p[i], p1);  
92.         t2 = multi(p2, polygon.p[i + 1], p1);  
93.         if(fabs(t1) < eps || fabs(t2) < eps){  
94.             if(fabs(t1) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i];      
95.             if(fabs(t2) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i + 1];  
96.         }  
97.         else if(t1 < 0 && t2 < 0){  
98.             new_polygon.p[new_polygon.n++] = polygon.p[i];    
99.             new_polygon.p[new_polygon.n++] = polygon.p[i + 1];  
100.         }  
101.         else if(t1 * t2  < 0){  
102.             l1 = lineFromSegment(p1, p2);  
103.             l2 = lineFromSegment(polygon.p[i], polygon.p[i + 1]);  
104.             interp = LineInter(l1, l2);   
105.             if(t1 < 0) {  
106.                 new_polygon.p[new_polygon.n++] = polygon.p[i];  
107.                 new_polygon.p[new_polygon.n++] = interp;  
108.             }  
109.             else {  
110.                 new_polygon.p[new_polygon.n++] = interp;  
111.                 new_polygon.p[new_polygon.n++] = polygon.p[i + 1];        
112.             }  
113.         }  
114.     }  
115.     polygon.n = 0;  
116.     if(new_polygon.n == 0) return polygon;  
117.     polygon.p[polygon.n++] = new_polygon.p[0];  
118.     for(i = 1;i < new_polygon.n;i++){  
119.         if(!same(new_polygon.p[i], new_polygon.p[i - 1])){  
120.             polygon.p[polygon.n++] = new_polygon.p[i];  
121.         }     
122.     }  
123.     if(polygon.n != 1 && same(polygon.p[polygon.n - 1], polygon.p[0])) polygon.n--;  
124.     polygon.p[polygon.n] = polygon.p[0];  
125.     return polygon;  
126. }  
127.   
128.   
129. double polygonArea(TPolygon p)  
130. {  
131.     //已知多边形各顶点的坐标，求其面积  
132.     int i, n;  
133.     double area;  
134.     n = p.n;  
135.     area = 0;  
136.     for(i = 1;i <= n;i++){  
137.         area += (p.p[i - 1].x * p.p[i % n].y - p.p[i % n].x * p.p[i - 1].y);  
138.     }    
139.     return area / 2;    
140. }  
141.   
142.   
143. void ChangeClockwise(TPolygon &polygon)  
144. {  
145.     TPoint tmp;  
146.     int i;  
147.     for(i = 0;i <= (polygon.n - 1) / 2;i++){  
148.         tmp = polygon.p[i];  
149.         polygon.p[i] = polygon.p[polygon.n - 1 - i];  
150.         polygon.p[polygon.n - 1 - i] = tmp;           
151.     }  
152. }  
153.   
154.   
155. int main()  
156. {  
157.     int test, i, j;  
158.     double area;  
159.     TPolygon polygon, new_polygon;  
160.     scanf("%d", &test);  
161.     while(test--){  
162.         scanf("%d", &polygon.n);  
163.         for(i = 0;i <= polygon.n - 1;i++){  
164.             scanf("%lf%lf", &polygon.p[i].x, &polygon.p[i].y);  
165.         }  
166.         /*若是逆时针转化为顺时针*/  
167.         if(polygonArea(polygon) > 0) ChangeClockwise(polygon);  
168.         polygon.p[polygon.n] = polygon.p[0];  
169.         new_polygon = polygon;  
170.         for(i = 0;i <= polygon.n - 1;i++){  
171.             new_polygon = Cut_polygon(polygon.p[i], polygon.p[i + 1], new_polygon);  
172.         }     
173.         area = polygonArea(new_polygon);  
174.         if(area < 0) printf("%.2lf\n", -area);  
175.         else printf("%.2lf\n", area);  
176.     }  
177.     return 0;     
178. }  
179.   
180.   
181. //是否存在  
182.   
183.   
184. #include <stdio.h>  
185. #include <math.h>  
186.   
187.   
188. #define Maxn 3005  
189. const double eps = 1e-10;  
190.   
191.   
192. typedef struct TPodouble   
193. {  
194.     double x;  
195.     double y;  
196. }TPoint;  
197.   
198.   
199. typedef struct TPolygon  
200. {  
201.     TPoint p[Maxn];  
202.     int n;  
203. };  
204.   
205.   
206. typedef struct TLine  
207. {  
208.     double a, b, c;  
209. }TLine;  
210.   
211.   
212. bool same(TPoint p1, TPoint p2)  
213. {  
214.     if(p1.x != p2.x) return false;  
215.     if(p1.y != p2.y) return false;  
216.     return true;  
217. }  
218.    
219. double multi(TPoint p1, TPoint p2, TPoint p0)  
220. {  
221.     //求矢量[p0, p1], [p0, p2]的叉积  
222.     //p0是顶点   
223.     return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);  
224.     //若结果等于0，则这三点共线  
225.     //若结果大于0，则p0p2在p0p1的逆时针方向  
226.     //若结果小于0，则p0p2在p0p1的顺时针方向   
227. }  
228.   
229.   
230. TLine lineFromSegment(TPoint p1, TPoint p2)  
231. {  
232.     //线段所在直线,返回直线方程的三个系统   
233.     TLine tmp;  
234.     tmp.a = p2.y - p1.y;  
235.     tmp.b = p1.x - p2.x;  
236.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
237.     return tmp;  
238. }  
239.   
240.   
241. TPoint LineInter(TLine l1, TLine l2)  
242. {  
243.     //求两直线得交点坐标  
244.     TPoint tmp;   
245.     double a1 = l1.a;  
246.     double b1 = l1.b;  
247.     double c1 = l1.c;  
248.     double a2 = l2.a;  
249.     double b2 = l2.b;  
250.     double c2 = l2.c;  
251.     //注意这里b1 = 0   
252.     if(fabs(b1) < eps){  
253.         tmp.x = -c1 / a1;    
254.         tmp.y = (-c2 - a2 * tmp.x) / b2;  
255.     }         
256.     else{  
257.         tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);  
258.         tmp.y = (-c1 - a1 * tmp.x) / b1;  
259.     }  
260.     return tmp;  
261. }  
262.   
263.   
264. TPolygon Cut_polygon(TPoint p1, TPoint p2, TPolygon polygon)  
265. {  
266.     TPolygon new_polygon;  
267.     TPoint interp;  
268.     TLine l1, l2;  
269.     int i, j;  
270.     double t1, t2;  
271.     new_polygon.n = 0;  
272.     for(i = 0;i <= polygon.n - 1;i++){  
273.         t1 = multi(p2, polygon.p[i], p1);  
274.         t2 = multi(p2, polygon.p[i + 1], p1);  
275.         if(fabs(t1) < eps || fabs(t2) < eps){  
276.             if(fabs(t1) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i];      
277.             if(fabs(t2) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i + 1];  
278.         }  
279.         else if(t1 < 0 && t2 < 0){  
280.             new_polygon.p[new_polygon.n++] = polygon.p[i];    
281.             new_polygon.p[new_polygon.n++] = polygon.p[i + 1];  
282.         }  
283.         else if(t1 * t2  < 0){  
284.             l1 = lineFromSegment(p1, p2);  
285.             l2 = lineFromSegment(polygon.p[i], polygon.p[i + 1]);  
286.             interp = LineInter(l1, l2);   
287.             if(t1 < 0) {  
288.                 new_polygon.p[new_polygon.n++] = polygon.p[i];  
289.                 new_polygon.p[new_polygon.n++] = interp;  
290.             }  
291.             else {  
292.                 new_polygon.p[new_polygon.n++] = interp;  
293.                 new_polygon.p[new_polygon.n++] = polygon.p[i + 1];        
294.             }  
295.         }  
296.     }  
297.     polygon.n = 0;  
298.     if(new_polygon.n == 0) return polygon;  
299.     polygon.p[polygon.n++] = new_polygon.p[0];  
300.     for(i = 1;i < new_polygon.n;i++){  
301.         if(!same(new_polygon.p[i], new_polygon.p[i - 1])){  
302.             polygon.p[polygon.n++] = new_polygon.p[i];  
303.         }     
304.     }  
305.     if(polygon.n != 1 && same(polygon.p[polygon.n - 1], polygon.p[0])) polygon.n--;  
306.     polygon.p[polygon.n] = polygon.p[0];  
307.     return polygon;  
308. }  
309.   
310.   
311. void ChangeClockwise(TPolygon &polygon)  
312. {  
313.     TPoint tmp;  
314.     int i;  
315.     for(i = 0;i <= (polygon.n - 1) / 2;i++){  
316.         tmp = polygon.p[i];  
317.         polygon.p[i] = polygon.p[polygon.n - 1 - i];  
318.         polygon.p[polygon.n - 1 - i] = tmp;           
319.     }  
320. }  
321.   
322.   
323. double polygonArea(TPolygon p)  
324. {  
325.     //已知多边形各顶点的坐标，求其面积  
326.     double area;  
327.     int i, n;  
328.     n = p.n;  
329.     area = 0;  
330.     for(i = 1;i <= n;i++){  
331.         area += (p.p[i - 1].x * p.p[i % n].y - p.p[i % n].x * p.p[i - 1].y);  
332.     }    
333.     return area / 2;     
334. }  
335.   
336.   
337. int main()  
338. {  
339.     int i, j;  
340.     TPolygon polygon, new_polygon;  
341.     while(scanf("%d", &polygon.n) && polygon.n){  
342.         for(i = 0;i <= polygon.n - 1;i++){  
343.             scanf("%lf%lf", &polygon.p[i].x, &polygon.p[i].y);  
344.         }  
345.         /*若是逆时针转化为顺时针*/  
346.         if(polygonArea(polygon) > 0) ChangeClockwise(polygon);         
347.         polygon.p[polygon.n] = polygon.p[0];  
348.         new_polygon = polygon;  
349.         for(i = 0;i <= polygon.n - 1;i++){  
350.             new_polygon = Cut_polygon(polygon.p[i], polygon.p[i + 1], new_polygon);  
351.         }  
352.         if(new_polygon.n > 0) printf("1\n");  
353.         else printf("0\n");   
354.     }  
355.     return 0;     
356. }  

1.16.8 pku_2600_二分+圆的参数方程

[cpp]  view plain  copy

1. #include <stdio.h>  
2. #include <math.h>  
3.   
4.   
5. const double eps = 1e-4;  
6. const double pi = acos(-1.0);  
7.   
8.   
9. struct TPoint   
10. {  
11.     double x, y;  
12. }p[60], a[60];  
13. double angle[60];  
14.   
15.   
16. double multi(TPoint p1, TPoint p2, TPoint p0)  
17. {  
18.     return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);  
19. }  
20.   
21.   
22. TPoint fine_a2(TPoint a1, TPoint m, double angle1)  
23. {  
24.     TPoint a2;  
25.     double r, angle2, angle3;  
26.     r = sqrt((a1.x - m.x) * (a1.x - m.x) + (a1.y - m.y) * (a1.y - m.y));  
27.     angle2 = acos((a1.x - m.x) / r);  
28.     if(a1.y < m.y) {  
29.         if(angle2 <= pi / 2) angle2 = -angle2;  
30.         if(angle2 > pi / 2) angle2 = 3 * pi / 2 - (angle2 - pi / 2);  
31.     }  
32.     angle3 = angle2 - angle1;  
33.     a2.x = m.x + r * cos(angle3);  
34.     a2.y = m.y + r * sin(angle3);  
35.     if(multi(m, a2, a1) < 0) return a2;  
36.     angle3 = angle2 + angle1;  
37.     a2.x = m.x + r * cos(angle3);  
38.     a2.y = m.y + r * sin(angle3);  
39.     if(multi(m, a2, a1) < 0) return a2;    
40. }  
41.   
42.   
43. int main()  
44. {  
45.     int n, i, j;  
46.     while(scanf("%d", &n) != EOF){  
47.         for(i = 0;i < n;i++){  
48.             scanf("%lf%lf", &p[i].x, &p[i].y);  
49.         }  
50.         for(i = 0;i < n;i++){  
51.             scanf("%lf", &angle[i]);  
52.             angle[i] = angle[i] * pi / 180;  
53.         }  
54.         a[0].x = 0;  
55.         a[0].y = 0;  
56.         while(1){  
57.             for(i = 1;i <= n;i++){  
58.                 a[i] = fine_a2(a[i - 1], p[i - 1], angle[i - 1]);  
59.             }  
60.             if(fabs(a[n].x - a[0].x) <= eps    
61.               && fabs(a[n].y - a[0].y) <= eps) break;  
62.             else {  
63.                 a[0].x = (a[0].x + a[n].x) / 2;  
64.                 a[0].y = (a[0].y + a[n].y) / 2;  
65.             }  
66.         }  
67.         for(i = 0;i < n;i++){  
68.             printf("%.0lf %.0lf\n", a[i].x, a[i].y);  
69.         }         
70.     }  
71.     return 0;  
72. }  
73. 1.16.9 pku_1151_矩形相交的面积  
74. /* 
75. 大牛的思想  
76. 题目给出 n 个矩形，要求它们的面积并。具体做法是离散化。 
77. 先把 2n 个 x 坐标排序去重，然后再把所有水平线段（ 
78. 要记录是矩形上边还是下边）按 y 坐标排序。 
79. 最后对于每一小段区间 (x[i], x[i + 1]) 扫描所有的水平线段， 
80. 求出这些水平线段在小区间内覆盖的面积。总的时间复杂度是 O(n^2)。 
81. 利用线段树，可以优化到 O(nlogn)。 
82. */  
83.   
84.   
85. #include <stdio.h>  
86. #include <math.h>  
87. #include <stdlib.h>  
88.   
89.   
90. #define up  1  
91. #define down -1  
92.   
93.   
94.   
95.   
96. typedef struct TSeg  
97. {  
98.     double l, r;  
99.     double y;  
100.     int UpOrDown;  
101. }TSeg;  
102. TSeg seg[210];  
103. int segn;  
104. double x[210];  
105. int xn;  
106.   
107.   
108. int cmp1(const void *a, const void *b)  
109. {  
110.     if(*(double *)a < *(double *)b) return -1;  
111.     else return 1;  
112. }  
113.   
114.   
115. int cmp2(const void *a, const void *b)  
116. {  
117.     TSeg *c = (TSeg *)a;  
118.     TSeg *d = (TSeg *)b;  
119.     if(c->y < d->y) return -1;  
120.     else return 1;  
121. }  
122.   
123.   
124. void movex(int t, int &xn)  
125. {  
126.     int i;  
127.     for(i = t;i <= xn - 1;i++){  
128.         x[i] = x[i + 1];  
129.     }  
130.     xn--;  
131. }  
132.   
133.   
134. int main()  
135. {  
136.     //freopen("in.in", "r", stdin);  
137.     //freopen("out.out", "w", stdout);  
138.     int n, i, j, cnt, test = 1;  
139.     double x1, y1, x2, y2, ylow, area;  
140.     while(scanf("%d", &n) != EOF && n){  
141.         xn = 0;  
142.         segn = 0;         
143.         for(i = 0;i < n;i++){  
144.             scanf("%lf%lf%lf%lf", &x1, &y1, &x2, &y2);  
145.             x[xn++] = x1;  
146.             x[xn++] = x2;  
147.             seg[segn].l = x1;  
148.             seg[segn].r = x2;         
149.             seg[segn].y = y1;  
150.             seg[segn++].UpOrDown = up;  
151.             seg[segn].l = x1;  
152.             seg[segn].r = x2;  
153.             seg[segn].y = y2;  
154.             seg[segn++].UpOrDown = down;              
155.         }  
156.         qsort(x, xn, sizeof(x[0]), cmp1);  
157.         /*除掉重复的x*/  
158.         for(i = 1;i < xn;){  
159.             if(x[i] == x[i - 1]) movex(i, xn);  
160.             else i++;  
161.         }  
162.         qsort(seg, segn, sizeof(seg[0]), cmp2);   
163.         area = 0.0;  
164.         for(i = 0;i < xn - 1;i++){  
165.             cnt = 0;          
166.             for(j = 0;j < segn;j++){  
167.                 if(seg[j].l <= x[i] && seg[j].r >= x[i + 1]){  
168.                     if(cnt == 0) ylow = seg[j].y;  
169.                     if(seg[j].UpOrDown == down) cnt++;  
170.                     else cnt--;  
171.                     if(cnt == 0) area += (x[i + 1] - x[i]) * (seg[j].y - ylow);                   
172.                 }  
173.             }  
174.         }  
175.         printf("Test case #%d\n", test++);  
176.         printf("Total explored area: %.2lf\n", area);         
177.     }  
178.     return 0;  
179. }  

1.16.10 pku_1118_共线最多的点的个数

[cpp]  view plain  copy

1. /* 
2. 2617120 chenhaifeng 1118 Accepted 512K 1890MS C++ 977B 2007-09-04 18:43:26  
3. 直接O(n^3)超时，用一个标记数组，标记i,j所做直线已经查找过，可以跳过 
4.  
5.  
6. 大牛的思想 
7. 朴素做法是 O(n3) 的，超时。我的做法是枚举每个点， 
8. 然后求其它点和它连线的斜率，再排序。这样就得到经过 
9. 该点的直线最多能经过几个点。求个最大值就行了。复 
10. 杂度是 O(n2logn) 的。把排序换成 hash， 
11. 可以优化到 O(n2)。  
12. 2617134 chenhaifeng 1118 Accepted 276K 312MS G++ 1394B 2007-09-04 18:49:08  
13. */  
14. #include <stdio.h>  
15. #include <math.h>  
16.   
17.   
18. bool f[705][705];  
19. int a[705];  
20.   
21.   
22. int main()  
23. {  
24.     int n, i, j, s, num, maxn;  
25.     int x[705], y[705];  
26.     int t, m;  
27.   
28.   
29.       
30.     while(scanf("%d", &n) != EOF && n){  
31.         for(i = 0;i <= n - 1;i++){  
32.             scanf("%d%d", &x[i], &y[i]);  
33.         }  
34.         maxn = -1;  
35.         for(i = 0;i <= n - 1;i++){  
36.             for(j = i;j <= n - 1;j++){  
37.                 f[i][j] = false;  
38.             }  
39.         }  
40.         for(i = 0;i <= n - 1;i++){  
41.             for(j = i + 1;j <= n - 1;j++){  
42.                 if(f[i][j] == true) continue;  
43.                 if(n - j < maxn) break;  
44.                 num = 2;  
45.                 t = 2;  
46.                 a[0] = i;  
47.                 a[1] = j;  
48.                 f[i][j] = true;   
49.                 for(s = j + 1;s <= n - 1;s++){  
50.                     if(f[i][s] == true || f[j][s] == true) continue;  
51.                     if((y[i] - y[s]) * (x[j] - x[s]) == (x[i] - x[s]) * (y[j] - y[s])){  
52.                          num++;   
53.                          a[t] = s;  
54.                          for(m = 0;m <= t - 1;m++){  
55.                                 f[m][s] = true;  
56.                         }  
57.                         t++;      
58.                     }                 
59.                 }  
60.                 if(num > maxn) maxn = num;     
61.             }             
62.         }  
63.         printf("%d\n", maxn);             
64.     }  
65.     return 0;  
66. }   

1.16.11 pku2826_线段围成的区域可储水量

[cpp]  view plain  copy

1. /* 
2. 两条线不相交， 
3. 左边或右边的口被遮住， 
4. 交点是某条线的那个纵坐标较高的那点 
5. 某条线段水平放置 
6. */  
7. #include <stdio.h>  
8. #include <math.h>  
9.   
10.   
11. #define eps 1e-8  
12.   
13.   
14. struct TPoint  
15. {  
16.     double x, y;  
17. };  
18. struct TLine  
19. {  
20.     double a, b, c;  
21. };  
22.   
23.   
24. int same(TPoint p1, TPoint p2)  
25. {  
26.     if(fabs(p1.x - p2.x) > eps) return 0;  
27.     if(fabs(p1.y - p2.y) > eps) return 0;  
28.     return 1;  
29. }  
30.   
31.   
32. double min(double x, double y)  
33. {  
34.     if(x < y) return x;  
35.     else return y;   
36. }  
37.   
38.   
39. double max(double x, double y)  
40. {  
41.     if(x > y) return x;  
42.     else return y;   
43. }  
44.   
45.   
46. double multi(TPoint p1, TPoint p2, TPoint p0)  
47. {   
48.     return (p1.x - p0.x) * (p2.y - p0.y)   
49.          - (p2.x - p0.x) * (p1.y - p0.y);  
50. }  
51.   
52.   
53. bool isIntersected(TPoint s1, TPoint e1, TPoint s2, TPoint e2)  
54. {  
55.     if(  
56.     (max(s1.x, e1.x) >= min(s2.x, e2.x)) &&  
57.     (max(s2.x, e2.x) >= min(s1.x, e1.x)) &&  
58.     (max(s1.y, e1.y) >= min(s2.y, e2.y)) &&  
59.     (max(s2.y, e2.y) >= min(s1.y, e1.y)) &&  
60.     (multi(s2, e1, s1) * multi(e1, e2, s1) >= 0) &&  
61.     (multi(s1, e2, s2) * multi(e2, e1, s2) >= 0)  
62.     )  return true;  
63.       
64.     return false;      
65. }  
66.   
67.   
68. TLine lineFromSegment(TPoint p1, TPoint p2)  
69. {  
70.     TLine tmp;  
71.     tmp.a = p2.y - p1.y;  
72.     tmp.b = p1.x - p2.x;  
73.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
74.     return tmp;  
75. }  
76.   
77.   
78. TPoint LineInter(TLine l1, TLine l2)  
79. {  
80.     TPoint tmp;   
81.     double a1 = l1.a;  
82.     double b1 = l1.b;  
83.     double c1 = l1.c;  
84.     double a2 = l2.a;  
85.     double b2 = l2.b;  
86.     double c2 = l2.c;  
87.     if(fabs(b1) < eps){  
88.         tmp.x = -c1 / a1;    
89.         tmp.y = (-c2 - a2 * tmp.x) / b2;  
90.     }         
91.     else{  
92.         tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);  
93.         tmp.y = (-c1 - a1 * tmp.x) / b1;  
94.     }  
95.     return tmp;  
96. }  
97.   
98.   
99. double triangleArea(TPoint p1, TPoint p2, TPoint p3)  
100. {  
101.     TPoint p4, p5;  
102.     p4.x = p2.x - p1.x;  
103.     p4.y = p2.y - p1.y;  
104.     p5.x = p3.x - p1.x;  
105.     p5.y = p3.y - p1.y;  
106.     return fabs(p5.x * p4.y - p5.y * p4.x) / 2;   
107. }  
108.   
109.   
110. double find_x(double y, TLine line)  
111. {  
112.     return (-line.c - line.b * y) / line.a;  
113. }  
114.   
115.   
116. double find_y(double x, TLine line)  
117. {  
118.     if(fabs(line.b) < eps)  
119.     {  
120.         return -1e250;  
121.     }  
122.     else   
123.     {  
124.         return (-line.c - line.a  * x) / line.b;  
125.     }  
126. }  
127.   
128.   
129. int main()  
130. {  
131.     //freopen("in.in", "r", stdin);  
132.     //freopen("out.out", "w", stdout);  
133.     int test;  
134.     double miny, y;  
135.     TLine l1, l2;  
136.     TPoint p1, p2, p3, p4, inter;  
137.     TPoint tp1, tp2;  
138.     scanf("%d", &test);  
139.     while(test--)  
140.     {  
141.         scanf("%lf%lf%lf%lf%lf%lf%lf%lf", &p1.x, &p1.y,   
142.         &p2.x, &p2.y, &p3.x, &p3.y, &p4.x, &p4.y);  
143.         if(same(p1, p2) || same(p3, p4)   
144.            || !isIntersected(p1, p2, p3, p4)  
145.            || fabs(p1.y - p2.y) < eps //平行与x轴   
146.            || fabs(p3.y - p4.y) < eps  
147.           )  
148.         {  
149.             printf("0.00\n");  
150.             continue;  
151.         }  
152.         l1 = lineFromSegment(p1, p2);  
153.         l2 = lineFromSegment(p3, p4);  
154.         inter = LineInter(l1, l2);  
155.         if(p1.y > p2.y) tp1 = p1;  
156.         else tp1 = p2;  
157.         if(p3.y > p4.y) tp2 = p3;  
158.         else tp2 = p4;  
159.         if(tp1.y < tp2.y)  
160.         {  
161.             if(tp1.x >= min(p4.x, p3.x) && tp1.x <= max(p4.x, p3.x))  
162.             {  
163.                 y = find_y(tp1.x, l2);  
164.                 if(y >= tp1.y)  
165.                 {  
166.                     printf("0.00\n");  
167.                     continue;  
168.                 }     
169.             }  
170.             miny = tp1.y;  
171.         }  
172.         else  
173.         {  
174.             if(tp2.x >= min(p1.x, p2.x) && tp2.x <= max(p1.x, p2.x))  
175.             {  
176.                 y = find_y(tp2.x, l1);  
177.                 if(y >= tp2.y)  
178.                 {  
179.                     printf("0.00\n");  
180.                     continue;  
181.                 }     
182.             }  
183.             miny = tp2.y;  
184.         }  
185.         if(fabs(miny - inter.y) < eps)  
186.         {  
187.             printf("0.00\n");  
188.             continue;             
189.         }  
190.         tp1.x = find_x(miny, l1);  
191.         tp2.x = find_x(miny, l2);  
192.         tp1.y = tp2.y = miny;  
193.         printf("%.2lf\n", triangleArea(tp1, tp2, inter));     
194.     }   
195.     return 0;  
196. }  

1.16.12 Pick公式

[cpp]  view plain  copy

1. //A = b / 2 + i -1  其中 b 与 i 分別表示在边界上及內部的格子点之个數  
2. //http://www.hwp.idv.tw/bbs1/htm/%A6V%B6q%B7L%BFn%A4%C0/%A6V%B6q%B7L%BFn%A4%C0.htm  
3. // http://acm.pku.edu.cn/JudgeOnline/problem?id=2954  
4.   
5.   
6. #include <stdio.h>  
7. #include <stdlib.h>  
8.   
9.   
10. typedef struct TPoint  
11. {  
12.     int x;  
13.     int y;  
14. }TPoint;  
15.   
16.   
17. typedef struct TLine  
18. {  
19.     int a, b, c;  
20. }TLine;  
21.   
22.   
23. int triangleArea(TPoint p1, TPoint p2, TPoint p3)  
24. {  
25.     //已知三角形三个顶点的坐标，求三角形的面积   
26.     int k = p1.x * p2.y + p2.x * p3.y + p3.x * p1.y  
27.       - p2.x * p1.y - p3.x * p2.y - p1.x * p3.y;    
28.     if(k < 0) return -k;  
29.     else return k;  
30. }  
31.   
32.   
33. TLine lineFromSegment(TPoint p1, TPoint p2)  
34. {  
35.     //线段所在直线,返回直线方程的三个系统   
36.     TLine tmp;  
37.     tmp.a = p2.y - p1.y;  
38.     tmp.b = p1.x - p2.x;  
39.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
40.     return tmp;  
41. }  
42.   
43.   
44. void swap(int &a, int &b)  
45. {  
46.     int t;  
47.     t = a;  
48.     a = b;  
49.     b = t;  
50. }  
51.   
52.   
53. int Count(TPoint p1, TPoint p2)  
54. {  
55.     int i, sum = 0, y;  
56.     TLine l1 =  lineFromSegment(p1, p2);  
57.     if(l1.b == 0) return abs(p2.y - p1.y) + 1;  
58.     if(p1.x > p2.x) swap(p1.x, p2.x); //这里没有交换WA两次   
59.     for(i = p1.x;i <= p2.x;i++){  
60.         y = -l1.c - l1.a * i;  
61.         if(y % l1.b == 0) sum++;  
62.     }  
63.     return sum;   
64. }  
65.   
66.   
67. int main()  
68. {  
69.     //freopen("in.in", "r", stdin);  
70.     //freopen("OUT.out", "w", stdout);  
71.     TPoint p1, p2, p3;  
72.     while(scanf("%d%d%d%d%d%d", &p1.x, &p1.y, &p2.x, &p2.y, &p3.x, &p3.y) != EOF){  
73.         if(p1.x == 0 && p1.y == 0 && p2.x == 0 && p2.y == 0 && p3.x == 0 && p3.y == 0) break;  
74.         int  A = triangleArea(p1, p2, p3);//A为面积的两倍   
75.         int b = 0;  
76.         int i;  
77.         b = Count(p1, p2) + Count(p1, p3) + Count(p3, p2) - 3;//3个顶点多各多加了一次   
78.         //i = A / 2- b / 2 + 1;  
79.         i = (A - b) / 2 + 1;  
80.         printf("%d\n", i);    
81.     }   
82.     return 0;     
83. }  

1.16.13 N点中三个点组成三角形面积最大

[cpp]  view plain  copy

1. //Rotating Calipers algorithm  
2.   
3.   
4.   
5.   
6. #include <stdio.h>  
7. #include <stdlib.h>  
8. #include <math.h>  
9.   
10.   
11. #define MaxNode 50005  
12.   
13.   
14. int stack[MaxNode];  
15. int top;  
16. double max;  
17.   
18.   
19. typedef struct TPoint  
20. {  
21.     int x;  
22.     int y;  
23. }TPoint;  
24. TPoint point[MaxNode];  
25.   
26.   
27. void swap(TPoint point[], int i, int j)  
28. {  
29.     TPoint tmp;  
30.     tmp = point[i];  
31.     point[i] = point[j];  
32.     point[j] = tmp;  
33. }  
34.   
35.   
36. double multi(TPoint p1, TPoint p2, TPoint p0)  
37. {  
38.     return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);  
39. }  
40.   
41.   
42. double distance(TPoint p1, TPoint p2)  
43. {  
44.     return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);  
45. }  
46.   
47.   
48. int cmp(const void *a, const void *b)  
49. {  
50.     TPoint *c = (TPoint *)a;  
51.     TPoint *d = (TPoint *)b;  
52.     double k = multi(*c, *d, point[0]);  
53.     if(k< 0) return 1;  
54.     else if(k == 0 && distance(*c, point[0]) >= distance(*d, point[0]))    
55.             return 1;  
56.     else return -1;     
57. }  
58.   
59.   
60. void grahamScan(int n)  
61. {   
62.     //Graham扫描求凸包  
63.     int i, u;   
64.       
65.     //将最左下的点调整到p[0]的位置  
66.     u = 0;  
67.     for(i = 1;i <= n - 1;i++){  
68.         if((point[i].y < point[u].y) ||   
69.              (point[i].y == point[u].y && point[i].x  < point[u].x))  
70.             u = i;        
71.     }   
72.     swap(point, 0, u);  
73.       
74.     //将平p[1]到p[n - 1]按按极角排序，可采用快速排序  
75.     qsort(point + 1, n - 1, sizeof(point[0]), cmp);  
76.       
77.     for(i = 0;i <= 2;i++) stack[i] = i;  
78.     top = 2;  
79.     for(i = 3;i <= n - 1;i++){  
80.         while(multi(point[i], point[stack[top]], point[stack[top - 1]]) >= 0){  
81.             top--;  
82.             if(top == 0) break;  
83.         }  
84.         top++;  
85.         stack[top] = i;  
86.     }  
87. }  
88.   
89.   
90. int main()  
91. {      
92.     double triangleArea(int i, int j, int k);  
93.     void PloygonTriangle();  
94.     int i, n;  
95.     while(scanf("%d", &n) && n != -1){  
96.         for(i = 0;i < n;i++)  
97.             scanf("%d%d", &point[i].x, &point[i].y);  
98.         if(n <= 2){  
99.             printf("0.00\n");  
100.             continue;        
101.         }  
102.         if(n == 3){  
103.             printf("%.2lf\n", triangleArea(0, 1, 2));  
104.             continue;   
105.         }  
106.         grahamScan(n);  
107.         PloygonTriangle();    
108.         printf("%.2lf\n", max);      
109.     }  
110.     return 0;  
111. }  
112.   
113.   
114. void PloygonTriangle()  
115. {  
116.     double triangleArea(int i, int j, int k);  
117.     int i, j , k;  
118.     double area, area1;  
119.     max = -1;  
120.     for(i = 0;i <= top - 2;i++){   
121.         k = -1;  
122.         for(j = i + 1; j <= top - 1;j++){   
123.             if(k <= j) k= j + 1;  
124.             area = triangleArea(stack[i], stack[j], stack[k]);  
125.             if(area > max) max = area;  
126.             while(k + 1 <= top){  
127.                 area1= triangleArea(stack[i], stack[j], stack[k + 1]);  
128.                 if(area1 < area) break;  
129.                 if(area1 > max) max = area1;  
130.                 area = area1;  
131.                 k++;  
132.             }  
133.         }  
134.     }  
135. }  
136.   
137.   
138. double triangleArea(int i, int j, int k)  
139. {  
140.     //已知三角形三个顶点的坐标，求三角形的面积   
141.     double l = fabs(point[i].x * point[j].y + point[j].x * point[k].y   
142.         + point[k].x * point[i].y - point[j].x * point[i].y   
143.         - point[k].x * point[j].y - point[i].x * point[k].y) / 2;    
144.     return l;  
145. }  

1.16.14 直线关于圆的反射

[cpp]  view plain  copy

1. /* 
2. fzu_1035 
3. 1.直线和圆的交点  
4. 2.点关于线的对称点 
5. 3.点到线的距离  
6. 4.直线方程  
7. */  
8. #include <iostream>  
9.   
10.   
11. #include <cmath>  
12.   
13.   
14. using namespace std;  
15.   
16.   
17. #define INF 999999999  
18. const double eps = 1e-6;  
19.   
20.   
21. int up;  
22.   
23.   
24. typedef struct TPoint  
25. {  
26.     double x;  
27.     double y;  
28. }TPoint;  
29.   
30.   
31. typedef struct TCircle  
32. {  
33.     TPoint center;  
34.     double r;  
35. }TCircle;  
36.   
37.   
38. typedef struct TLine  
39. {  
40.     //直线标准式中的系数   
41.     double a, b, c;  
42. }TLine;  
43.   
44.   
45. void SloveLine(TLine &line, TPoint start, TPoint dir)  
46. {  
47.     //根据直线上一点和直线的方向求直线的方程   
48.     if(dir.x == 0){  
49.         line.a = 1;  
50.         line.b = 0;  
51.         line.c = start.x;  
52.     }  
53.     else {  
54.         double k = dir.y / dir.x;  
55.         line.a = k;  
56.         line.b = -1;  
57.         line.c = start.y - k * start.x;          
58.     }      
59. }  
60.   
61.   
62. TLine lineFromSegment(TPoint p1, TPoint p2)  
63. {  
64.     //线段所在直线,返回直线方程的三个系统   
65.     TLine tmp;  
66.     tmp.a = p2.y - p1.y;  
67.     tmp.b = p1.x - p2.x;  
68.     tmp.c = p2.x * p1.y - p1.x * p2.y;  
69.     return tmp;  
70. }  
71.   
72.   
73. TPoint symmetricalPointofLine(TPoint p, TLine L)  
74. {  
75.     //p点关于直线L的对称点   
76.     TPoint p2;  
77.     double d;  
78.     d = L.a * L.a + L.b * L.b;  
79.     p2.x = (L.b * L.b * p.x - L.a * L.a * p.x -   
80.             2 * L.a * L.b * p.y - 2 * L.a * L.c) / d;  
81.     p2.y = (L.a * L.a * p.y - L.b * L.b * p.y -   
82.             2 * L.a * L.b * p.x - 2 * L.b * L.c) / d;  
83.     return p2;  
84. }  
85.   
86.   
87. double distanc(TPoint p1, TPoint p2)  
88. {  
89.     //计算平面上两个点之间的距离  
90.    return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));      
91. }  
92.   
93.   
94. bool samedir(TPoint dir, TPoint start, TPoint point)  
95. {  
96.     //判断方向   
97.     TPoint tmp;  
98.     tmp.x = point.x - start.x;  
99.     tmp.y = point.y - start.y;  
100.     if(tmp.x != 0 && dir.x != 0){  
101.         if(tmp.x / dir.x > 0) return true;  
102.         else return false;  
103.     }  
104.     else if(tmp.y != 0 && dir.y != 0){  
105.         if(tmp.y / dir.y > 0) return true;  
106.         else return false;   
107.     }  
108.     return true;   
109. }  
110.   
111.   
112. bool Intersected(TPoint &point, TLine line, const TCircle circle[],   
113.                  TPoint start, TPoint dir, int which)  
114. {  
115.     //如果圆与直线有(有效交点)交点就存放在变量point中   
116.     double a = line.a, b = line.b, c = line.c;  
117.     double x0 = circle[which].center.x, y0 = circle[which].center.y;  
118.     double r =  circle[which].r;  
119.     //有交点，求交点   
120.     double x2front = b * b + a * a;  
121.     double x1front = -2 * x0 * b * b + 2 * a * b * y0 + 2 * a * c;  
122.     double front = x0 * x0 * b * b + y0 * y0 * b * b   
123.         + c * c + 2 * c * y0 * b - b * b * r * r;  
124.     double d = x1front * x1front - 4 * x2front * front;  
125.     TPoint p1, p2;  
126.     bool k1, k2;  
127.     if(fabs(d) < eps){  
128.         //x2front不可能等于零   
129.        point.x = -x1front / x2front / 2;  
130.        point.y = (-c - a * point.x) / b;  
131.        //判断方向  
132.        if(samedir(dir, start, point)) return true;  
133.        else return false;  
134.     }  
135.     else if(d < 0) return false;  
136.     else {  
137.         p1.x = (-x1front + sqrt(d)) / 2 / x2front;  
138.         p1.y = (-c - a * p1.x) / b;  
139.         p2.x = (-x1front - sqrt(d)) / 2 / x2front;  
140.         p2.y = (-c - a * p2.x) / b;  
141.         k1 = samedir(dir, start, p1);  
142.         k2 = samedir(dir, start, p2);  
143.         if(k1 == false && k2 == false) return false;  
144.         if(k1 == true && k2 == true){  
145.             double dis1 = distanc(p1, start);  
146.             double dis2 = distanc(p2, start);  
147.             if(dis1 < dis2) point = p1;  
148.             else point = p2;  
149.             return true;  
150.         }  
151.         else if(k1 == true) point = p1;  
152.         else point = p2;  
153.         return true;  
154.     }     
155. }  
156.   
157.   
158. void Reflect(int &num, TCircle circle[], TPoint start, TPoint dir, int n)  
159. {  
160.     //反复反射   
161.     int i;  
162.     TLine line;  
163.     TPoint interpoint, newstart;  
164.     int u;  
165.     SloveLine(line, start, dir);  
166.     int tag = 0;  
167.     double mindis = INF;  
168.     for(i = 1;i <= n;i++){  
169.         if(i != up && Intersected(interpoint, line, circle, start, dir, i)){  
170.             double dis = distanc(start, interpoint);  
171.             if(dis < mindis){  
172.                 tag = 1;  
173.                 u = i;    
174.                 mindis = dis;  
175.                 newstart = interpoint;  
176.             }              
177.         }  
178.     }  
179.     if(tag == 0){  
180.         cout << "inf" << endl;  
181.         return ;  
182.     }  
183.     else {  
184.         if(num == 10){  
185.             cout << "..." << endl;  
186.            return ;  
187.         }  
188.         cout << u << " ";   
189.         num++;  
190.         //新的方向   
191.         TLine line1;  
192.         TPoint p;  
193.         line1 = lineFromSegment(newstart, circle[u].center);  
194.         if(fabs(line1.a * start.x + line1.b * start.y +line1.c) <= eps){  
195.             dir.x = -dir.x;  
196.             dir.y = -dir.y;  
197.         }  
198.         else {  
199.             p = symmetricalPointofLine(start, line1);//start的对称点   
200.             dir.x = p.x - newstart.x;  
201.             dir.y = p.y - newstart.y;   
202.         }  
203.           
204.         start = newstart;       
205.         up = u;  
206.         Reflect(num, circle, start, dir, n);       
207.     }   
208. }  
209.   
210.   
211. int main()  
212. {  
213.     //freopen("fzu_1035.in", "r", stdin);  
214.     //freopen("fzu_1035.out", "w", stdout);  
215.     int n, i, j, num, test = 1;  
216.     TCircle circle[30];  
217.     TPoint start, dir;  
218.     while(cin >> n && n){  
219.         for(i = 1;i <= n;i++){  
220.             cin >> circle[i].center.x >> circle[i].center.y >> circle[i].r;  
221.         }  
222.         cin >> start.x >> start.y >> dir.x >> dir.y;  
223.           
224.         cout << "Scene " << test++ << endl;  
225.   
226.   
227.         num = 0;  
228.         up = -1;  
229.         Reflect(num, circle, start, dir, n);  
230.         cout << endl;  
231.     }   
232.     return 0;  
233. }  

1.16.15 pku2002_3432_N个点最多组成多少个正方形(hao)

[cpp]  view plain  copy

1. #include <stdio.h>   
2. #include <stdlib.h>   
3. #include <math.h>   
4. #define eps 1e-6  
5. #define pi acos(-1.0)  
6.   
7.   
8. #define PRIME 9991  
9.   
10.   
11. struct point  
12. {   
13.     int x, y;   
14. }p[2201];   
15. int n;  
16.   
17.   
18. struct HASH  
19. {  
20.     int cnt;  
21.     int next;  
22. }hash[50000];  
23. int hashl;  
24.   
25.   
26. int Hash(int n)  
27. {  
28.     int i = n % PRIME;  
29.     while(hash[i].next != -1){  
30.         if(hash[hash[i].next].cnt == n) return 1;  
31.         else if(hash[hash[i].next].cnt > n) break;  
32.         i = hash[i].next;  
33.     }  
34.     hash[hashl].cnt = n;  
35.     hash[hashl].next = hash[i].next;  
36.     hash[i].next = hashl;  
37.     hashl++;   
38.     return 0;  
39. }  
40.   
41.   
42. int Hash2(int n)  
43. {  
44.     int i = n % PRIME;  
45.     while(hash[i].next != -1){  
46.         if(hash[hash[i].next].cnt == n) return 1;  
47.         else if(hash[hash[i].next].cnt > n) return 0;  
48.         i = hash[i].next;  
49.     }  
50.     return 0;  
51. }  
52.   
53.   
54. int check(double ax, double ay, int &x, int &y)  
55. {  
56.     int a0 = (int)ax;  
57.     int b0 = (int)ay;  
58.     int tag1 = 0, tag2 = 0;     
59.     if(fabs(a0 - ax) < eps){  
60.         tag1 = 1;  
61.         x = a0;  
62.     }  
63.     else if(fabs(a0 + 1 - ax) < eps){  
64.         tag1 = 1;  
65.         x = a0 + 1;  
66.     }  
67.     if(fabs(b0 - ay) < eps){  
68.         tag2 = 1;  
69.         y = b0;  
70.     }  
71.     else if(fabs(b0 + 1 - ay) < eps){  
72.         y = b0 + 1;  
73.         tag2 = 1;  
74.     }  
75.     if(tag1 == 1 && tag2 == 1) return 1;  
76.     else return 0;  
77. }  
78.   
79.   
80. int squares(point p1, point p2, point &p3, point &p4)  
81. {  
82.     double a = (double)p2.x - p1.x;  
83.     double b = (double)p2.y - p1.y;  
84.     double midx = ((double)p1.x + p2.x) / 2;  
85.     double midy = ((double)p1.y + p2.y) / 2;  
86.     double tmp = a * a + b * b;  
87.     double x1 = sqrt(b * b) / 2;  
88.     double y1;  
89.     if(fabs(b) < eps) y1 = sqrt(a * a + b * b) / 2;  
90.     else y1 = -a * x1 / b;  
91.     x1 += midx;  
92.     y1 += midy;  
93.     if(check(x1, y1, p3.x, p3.y) == 0) return 0;      
94.     x1 = 2 * midx - x1;  
95.     y1 = 2 * midy - y1;   
96.     if(check(x1, y1, p4.x, p4.y) == 0) return 0;   
97.     return 1;  
98. }  
99.   
100.   
101. int main()   
102. {   
103.     int i, j, cnt;   
104.     while(scanf("%d", &n) != EOF && n)  
105.     {  
106.         for(i = 0;i < PRIME;i++) hash[i].next = -1;  
107.         hashl = PRIME;  
108.         int x1, y1, x2, y2;  
109.         for (i = 0; i < n; i++){  
110.             scanf("%d%d", &p[i].x, &p[i].y);   
111.             Hash((p[i].x + 100000) * 100000 + p[i].y + 100000);  
112.         }     
113.         cnt = 0;   
114.         for (i = 0; i < n; i++){   
115.             for (j = i + 1; j < n; j++)   
116.             {   
117.                 point a, b;   
118.                 if(squares(p[i], p[j], a, b) == 0) continue;   
119.                 if(Hash2((a.x + 100000) * 100000 + a.y + 100000) == 0) continue;  
120.                 if(Hash2((b.x + 100000) * 100000 + b.y + 100000) == 0) continue;  
121.                 cnt++;   
122.             }   
123.         }  
124.         printf("%d\n", cnt / 2);   
125.     }  
126.     return 0;   
127. }   

1.16.16 pku1981_单位圆覆盖最多点(poj1981)CircleandPoints

[cpp]  view plain  copy

1. /* 
2. 平面上N个点，用一个半径R的圆去覆盖，最多能覆盖多少个点？ 
3.  
4.  
5. 比较经典的题目。 
6.  
7.  
8. 对每个点以R为半径画圆，对N个圆两两求交。这一步O(N^2)。问题转化为求被覆盖次数最多的弧。 
9.  
10.  
11. 对每一个圆，求其上的每段弧重叠次数。假如A圆与B圆相交。A上[PI/3, PI/2]的区间被B覆盖(PI为圆周率)。那么对于A圆，我们在PI/3处做一个+1标记，在PI/2处做一个-1标记。 
12.  
13.  
14. 对于[PI*5/3, PI*7/3]这样横跨0点的区间只要在0点处拆成两段即可。 
15.  
16.  
17. 将一个圆上的所有标记排序，从头开始扫描。初始ans = 0，碰到+1标记给ans++，碰到-1标记ans--。扫描过程中ans的最大值就是圆上被覆盖最多的弧。求所有圆的ans的最大值就是答案。 
18.  
19.  
20. 总复杂度O(N^2 * logN) 
21.  
22.  
23. */#include <stdio.h>  
24. #include <math.h>  
25.   
26.   
27. #define eps 1e-6  
28.   
29.   
30. struct point   
31. {  
32.     double x, y;  
33. };  
34.   
35.   
36. double dis(point p1, point p2)  
37. {  
38.     point p3;  
39.     p3.x = p2.x - p1.x;  
40.     p3.y = p2.y - p1.y;  
41.     return p3.x * p3.x + p3.y * p3.y;  
42. }  
43.   
44.   
45. point find_centre(point p1, point p2)  
46. {  
47.     point p3, mid, centre;  
48.     double b, c, ang;  
49.     p3.x = p2.x - p1.x;  
50.     p3.y = p2.y - p1.y;  
51.     mid.x = (p1.x + p2.x) / 2;  
52.     mid.y = (p1.y + p2.y) / 2;        
53.     b = dis(p1, mid);  
54.     c = sqrt(1 - b);  
55.     if(fabs(p3.y) < eps)//垂线的斜角90度  
56.     {  
57.         centre.x = mid.x;  
58.         centre.y = mid.y + c;  
59.     }  
60.     else   
61.     {  
62.         ang = atan(-p3.x / p3.y);  
63.         centre.x = mid.x + c * cos(ang);  
64.         centre.y = mid.y + c * sin(ang);  
65.     }  
66.     return centre;    
67. }  
68.   
69.   
70. int main()  
71. {  
72.     int n, ans, tmpans, i, j, k;  
73.     point p[305], centre;  
74.     double tmp;  
75.     while(scanf("%d", &n) && n)  
76.     {  
77.         for(i = 0;i < n;i++)  
78.             scanf("%lf%lf", &p[i].x, &p[i].y);  
79.         ans = 1;  
80.         for(i = 0;i < n;i++)  
81.             for(j = i + 1;j < n;j++)  
82.             {  
83.                 if(dis(p[i], p[j]) > 4) continue;  
84.                 tmpans = 0;  
85.                 centre = find_centre(p[i], p[j]);  
86.                 for(k = 0;k < n;k++)  
87.                 {  
88.                     //if(tmpans + n - k <= ans) break;  
89.                     tmp = dis(centre, p[k]);  
90.                     //if(tmp < 1.0 || fabs(tmp - 1.0) < eps) tmpans++;  
91.                     if(tmp <= 1.000001) tmpans++;  
92.                 }  
93.                 if(ans < tmpans) ans = tmpans;  
94.             }  
95.         printf("%d\n", ans);          
96.     }  
97.     return 0;  
98. }  

1.16.17 pku3668_GameofLine_N个点最多确定多少互不平行的直线(poj3668)

[cpp]  view plain  copy

1. #include <math.h>  
2. #include <stdio.h>  
3. #include <stdlib.h>  
4.   
5.   
6. #define eps 1e-6  
7. #define pi acos(-1)  
8.   
9.   
10. struct point   
11. {  
12.     double x, y;  
13. };  
14.   
15.   
16. double FindSlewRate(point p1, point p2)  
17. {  
18.     point p;  
19.     p.x = p2.x - p1.x;  
20.     p.y = p2.y - p1.y;  
21.     if(fabs(p.x) < eps) return pi / 2;  
22.     double tmp = atan(p.y / p.x);  
23.     if(tmp < 0) return pi + tmp;  
24.     return tmp;  
25. }  
26.   
27.   
28. int cmp(const void *a, const void *b)  
29. {  
30.     double *c = (double *)a;  
31.     double *d = (double *)b;  
32.     if(*c < *d) return -1;  
33.     return 1;  
34. }  
35.   
36.   
37. int main()  
38. {  
39.     int n, rt;  
40.     point p[205];  
41.     double rate[40005];  
42.     while(scanf("%d", &n) != EOF)  
43.     {  
44.         for(int i = 0;i < n;i++)  
45.             scanf("%lf%lf", &p[i].x ,&p[i].y);  
46.         rt = 0;  
47.         for(int i = 0;i < n;i++)  
48.             for(int j = i + 1;j < n;j++)  
49.                 rate[rt++] = FindSlewRate(p[i], p[j]);  
50.         qsort(rate, rt, sizeof(rate[0]), cmp);    
51.         int ans = 1;  
52.         for(int i = 1;i < rt;i++)  
53.             if(rate[i] > rate[i - 1]) ans++;  
54.         //注意这里写fabs(rate[i] - rate[i - 1]) > eps Wrong Answer   
55.         printf("%d\n", ans);  
56.     }   
57.     return 0;  
58. }  

1.16.18 求凸多边形直径

[cpp]  view plain  copy

1. #include<stdio.h>   
2. #include<math.h>   
3.   
4.   
5. #define eps 1e-6   
6. #define MaX 6000   
7.   
8.   
9. /*-----------多边形结构------------*/   
10. struct POLYGON{   
11.    int n;                                  //多边形顶点数   
12.    double x[MaX],y[MaX];        //顶点坐标   
13. }poly;   
14.   
15.   
16. int zd[100000][2],znum;      //跖对的集合和跖对的数量   
17.   
18.   
19.   
20.   
21. /*------------辅助函数-----------*/   
22. double dist(int a,int b,int c)   
23. {   
24.    double vx1,vx2,vy1,vy2;   
25.    vx1=poly.x[b]-poly.x[a]; vy1=poly.y[b]-poly.y[a];   
26.    vx2=poly.x[c]-poly.x[a]; vy2=poly.y[c]-poly.y[a];   
27.    return fabs(vx1*vy2 - vy1*vx2);   
28. }   
29.   
30.   
31.   
32.   
33. /*-------------求凸多边形直径的函数-------------*/   
34. double DIAMETER()   
35. {   
36.    znum=0;   
37.    int i,j,k=1;   
38.    double m,tmp;   
39.    while(dist(poly.n-1,0,k+1) > dist(poly.n-1,0,k)+eps)   
40.       k++;   
41.    i=0; j=k;   
42.    while(i<=k && j<poly.n)   
43.    {   
44.       zd[znum][0]=i; zd[znum++][1]=j;   
45.       while(dist(i,i+1,j+1)>dist(i,i+1,j)-eps && j<poly.n-1)   
46.       {   
47.          zd[znum][0]=i; zd[znum++][1]=j;   
48.          j++;   
49.       }   
50.       i++;   
51.    }   
52.    m=-1;   
53.    for(i=0;i<znum;i++)   
54.    {   
55.       tmp =(poly.x[zd[i][0]]-poly.x[zd[i][1]]) * (poly.x[zd[i][0]]-poly.x[zd[i][1]]);   
56.       tmp+=(poly.y[zd[i][0]]-poly.y[zd[i][1]]) * (poly.y[zd[i][0]]-poly.y[zd[i][1]]);   
57.       if(m<tmp) m=tmp;   
58.    }   
59.    return sqrt(m);   
60. }   
61.   
62.   
63. /*----------主函数----------*/   
64. int main()   
65. {   
66.    int i;   
67.    while(scanf("%d",&poly.n)==1)   
68.    {   
69.       for(i=0;i<poly.n;i++)   
70.          scanf("%lf %lf",&poly.x[i],&poly.y[i]);   
71.       printf("%.3lf\n",DIAMETER());   
72.    }   
73.    return 0;   
74. }  
75. 1.16.19 矩形面积并，周长并  
76. 见附录一5.11,5.12  
77. 1.16.20 pku2069 最小球覆盖   
78. 见附录一5.13,5.14  
79. //最小闭合球  
80. #include<stdio.h>  
81. #include<math.h>  
82. #include<memory>  
83. #include<stdlib.h>  
84. using namespace std;  
85. const double eps = 1e-10;  
86. struct point_type { double x, y, z; };  
87. int npoint, nouter ;  
88. point_type point [1000], outer[4], res;  
89. double radius, tmp ;  
90.   
91.   
92. inline double dist(point_type p1 , point_type p2)  
93. {  
94. double dx=p1.x-p2.x, dy=p1.y-p2.y,dz=p1.z-p2.z ;  
95. return ( dx*dx + dy*dy + dz*dz ) ;  
96. }  
97.   
98.   
99. inline double dot( point_type p1 , point_type p2 )  
100. {  
101. return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;  
102. }  
103.   
104.   
105. void ball()  
106. {  
107.     point_type q[3];   
108.     double m[3][3],sol[3],L[3],det; int i,j;  
109.     res.x=res.y=res.z=-1000;  
110.     radius=0;  
111.     switch ( nouter )   
112.         {  
113.         case 1 : res=outer[0]; break;  
114.         case 2 :  
115.                  res.x=(outer[0].x+outer[1].x)/2;  
116.                  res.y=(outer[0].y+outer[1].y)/2;  
117.                  res.z=(outer[0].z+outer[1].z)/2;  
118.                  radius=dist(res,outer[0]);  
119.                  break;  
120.         case 3 :  
121.                  for ( i=0; i<2; ++i ) {  
122.                     q[i].x=outer[i+1].x-outer[0].x;  
123.                     q[i].y=outer[i+1].y-outer[0].y;  
124.                     q[i].z=outer[i+1].z-outer[0].z;  
125.                     }  
126.                  for ( i=0; i<2; ++i )  
127.                      for ( j=0; j<2; ++j )  
128.                         m[i][j]=dot(q[i],q[j])*2 ;  
129.                  for ( i=0; i<2; ++i ) sol[i]=dot(q[i],q[i]);  
130.                  if (fabs(det=m[0][0]*m[1][1]-m[0][1]*m[1][0]) < eps ) return ;  
131.                    
132.                  L[0]=(sol[0]*m[1][1]-sol[1]*m[0][1])/det;  
133.                  L[1]=(sol[1]*m[0][0]-sol[0]*m[1][0])/det;  
134.                  res.x=outer[0].x+q[0].x*L[0]+q[1].x*L[1];  
135.                  res.y=outer[0].y+q[0].y*L[0]+q[1].y*L[1];  
136.                  res.z=outer[0].z+q[0].z*L[0]+q[1].z*L[1];  
137.                  radius=dist(res,outer[0]);  
138.                  break;  
139.         case 4 :  
140.                  for ( i=0; i<3; ++i ){  
141.                     q[i].x=outer[i+1].x-outer[0].x;  
142.                     q[i].y=outer[i+1].y-outer[0].y;  
143.                     q[i].z=outer[i+1].z-outer[0].z;  
144.                     sol[i]=dot(q[i],q[i]);  
145.                     }  
146.                  for ( i=0; i<3; ++i)   
147.                      for ( j=0; j<3; ++j) m[i][j]=dot(q[i],q[j])*2;  
148.                  det= m[0][0]*m[1][1]*m[2][2] + m[0][1]*m[1][2]*m[2][0]  
149.                      +m[0][2]*m[2][1]*m[1][0] - m[0][2]*m[1][1]*m[2][0]  
150.                      -m[0][1]*m[1][0]*m[2][2] - m[0][0]*m[1][2]*m[2][1];  
151.   
152.   
153.                  if ( fabs( det )<eps ) return;  
154.   
155.   
156.                  for ( j=0; j<3; ++j ){  
157.                     for ( i=0; i<3; ++i ) m[i][j]=sol[i];  
158.                      L[j]=( m[0][0]*m[1][1]*m[2][2] + m[0][1]*m[1][2]*m[2][0]  
159.                           + m[0][2]*m[2][1]*m[1][0] - m[0][2]*m[1][1]*m[2][0]  
160.                           - m[0][1]*m[1][0]*m[2][2] - m[0][0]*m[1][2]*m[2][1]  
161.                           ) / det;  
162.                     for( i=0; i<3; ++i ) m[i][j]=dot(q[i],q[j])*2;  
163.                     }  
164.                  res=outer[0];  
165.                  for ( i=0; i<3; ++i ) {  
166.                     res.x+=q[i].x*L[i];  
167.                     res.y+=q[i].y*L[i];  
168.                     res.z+=q[i].z*L[i];  
169.                     }  
170.                  radius=dist(res,outer[0]);  
171.         }  
172. }  
173.   
174.   
175. void minball(int n)  
176. {  
177.   ball();  
178.   if ( nouter <4 )  
179.        for ( int i=0; i<n; ++i )  
180.            if( dist(res,point[i])-radius>eps)  
181.             {  
182.                 outer[nouter]=point[i];  
183.                 ++nouter;  
184.                 minball(i);  
185.                 --nouter;  
186.                 if(i>0)  
187.                 {  
188.                   point_type Tt = point[i] ;  
189.                   memmove(&point[1], &point[0] , sizeof ( point_type )*i );  
190.                   point[0]=Tt;  
191.                 }  
192.             }  
193. }  
194.   
195.   
196.   
197.   
198.   
199.   
200. int main()  
201. {  
202.  int i;  
203.  while(scanf("%d",&npoint)!=EOF,npoint)  
204.     {  
205.      for(i=0;i<npoint;i++)  
206.        scanf("%lf%lf%lf",&point[i].x,&point[i].y,&point[i].z);  
207.      nouter=0;   
208.      minball(npoint);  
209.      printf("%.8lf\n",sqrt(radius)+eps);  
210.     }  
211.  return 0;  
212. }  

1.16.21 最大空凸包、最大空矩形

见附录一5.7, 6.15

1.16.22 求圆和多边形的交

[cpp]  view plain  copy

1. /*圆和简单多边形*/  
2. #include <cstdio>  
3. #include <cstring>  
4. #include <algorithm>  
5. #include <cmath>  
6. #include <cstdlib>  
7. #include <iostream>  
8. #include <ctime>  
9. using namespace std;  
10. #define M 30  
11. #define eps 1e-7  
12. const double PI = acos(-1.0);  
13.   
14.   
15. class pnt_type  
16. {  
17. public:  
18.     double x,y;  
19. };  
20. class state_type  
21. {  
22. public:  
23.     double angle;  
24.     double CoverArea;  
25. };  
26.   
27.   
28. pnt_type pnt[M];  
29. pnt_type center;  
30.   
31.   
32. int n;  
33. double R;  
34.   
35.   
36. bool read_data()  
37. {  
38.     n = 3;  
39.     int i;  
40.     if (cin >> pnt[1].x >> pnt[1].y)  
41.     {  
42.         for (i=2;i<=n;i++) cin >> pnt[i].x >> pnt[i].y;  
43.         cin >> center.x >> center.y >> R;  
44.         return true;  
45.     }  
46.     return false;  
47. }  
48. inline double Area2(pnt_type &a,pnt_type &b,pnt_type &c)  
49. {  
50.     return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);  
51. }  
52. inline double dot(pnt_type &a,pnt_type &b,pnt_type &c)  
53. {  
54.     return (b.x - a.x) * (c.x - a.x) + (b.y - a.y) * (c.y - a.y);  
55. }  
56. inline double dist(pnt_type &a,pnt_type &b)  
57. {  
58.     return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));  
59. }  
60.   
61.   
62. void init()  
63. {  
64.     int i;  
65.     double temp,sum;  
66.     for (i=2;i<n;i++)  
67.     {  
68.         temp = Area2(pnt[1],pnt[i],pnt[i + 1]);  
69.         sum += temp;  
70.     }  
71.     if (sum < 0) reverse(pnt + 1,pnt + n + 1);  
72.     pnt[n + 1] = pnt[1];  
73. }  
74.   
75.   
76. inline bool inCircle(pnt_type &s)  
77. {  
78.     return dist(center,s) <= R;  
79. }  
80.   
81.   
82. bool SameSide(pnt_type a,pnt_type b)  
83. {  
84.     if (dist(a,center) > dist(b,center)) swap(a,b);  
85.     return dot(a,b,center) < eps;  
86. }  
87.   
88.   
89. double ShadomOnCircle(pnt_type a,pnt_type b)  
90. {  
91.     double flag = Area2(center,a,b),res = 0;  
92.     if (fabs(flag) < eps) return 0;  
93.   
94.   
95.     bool ina = inCircle(a),inb = inCircle(b);  
96.     if (ina && inb)  
97.     {  
98.         res = fabs(Area2(center,a,b)) / 2;  
99.     }  
100.     else if (!ina && !inb)  
101.     {  
102.         if (SameSide(a,b))  
103.         {  
104.             double theta = acos(dot(center,a,b) / dist(center,a) / dist(center,b));  
105.             res = R * R * theta / 2;  
106.         }  
107.         else  
108.         {  
109.             double height = fabs(Area2(center,a,b)) / dist(a,b);  
110.             double theta = acos(dot(center,a,b) / dist(center,a) / dist(center,b));  
111.             if (height >= R)  
112.             {  
113.                 res = R * R * theta / 2;  
114.             }  
115.             else  
116.             {  
117.                 double _theta = 2 * acos(height / R);  
118.                 res = R * R * (theta - _theta) / 2 + R * R * sin(_theta) / 2;  
119.             }  
120.         }  
121.     }  
122.     else  
123.     {  
124.         if (!ina && inb) swap(a,b);  
125.         double height = fabs(Area2(center,a,b)) / dist(a,b);  
126.         double temp = dot(a,center,b);  
127.         double theta = acos(dot(center,a,b) / dist(center,a) / dist(center,b)),theta1,theta2;  
128.         if (fabs(temp) < eps)  
129.         {  
130.             double _theta = acos(height / R);  
131.             res += R * height / 2 * sin(_theta);  
132.             res += R * R / 2 * (theta - _theta);  
133.         }  
134.         else  
135.         {  
136.             theta1 = asin(height / R); theta2 = asin(height / dist(a,center));  
137.             if (temp > 0)  
138.             {  
139.                 res += dist(center,a) * R / 2 * sin(PI - theta1 - theta2);  
140.                 res += R * R / 2 * (theta + theta1 + theta2 - PI);  
141.             }  
142.             else  
143.             {  
144.                 res += dist(center,a) * R / 2 * sin(theta2 - theta1);  
145.                 res += R * R / 2 * (theta - theta2 + theta1);  
146.             }  
147.         }  
148.     }  
149.     if (flag < 0) return -res; else return res;  
150. }  
151.   
152.   
153. double Cover()  
154. {  
155.     int i;  
156.     double res = 0;  
157.     for (i=1;i<=n;i++)   
158.         res += ShadomOnCircle(pnt[i],pnt[i + 1]);  
159.     return res;  
160. }  
161.   
162.   
163. int main()  
164. {  
165.     double ans;  
166.     while (read_data())  
167.     {  
168.         init();  
169.         ans = Cover();  
170.         printf("%.2lf\n",ans);  
171.     }  
172.     return 0;  
173. }  

半平面交

[cpp]  view plain  copy

1. //Nlgn  
2. #include <stdio.h>  
3. #include <string.h>  
4. #include <math.h>  
5. #include <algorithm>  
6. using namespace std;  
7. #define maxn 20005  
8. #define eps 1e-10  
9.   
10.   
11. struct point  
12. {double x,y;};  
13. struct line  
14. {point s,e;double k;};  
15.   
16.   
17. line L[maxn];  
18. point S[maxn];  
19.   
20.   
21. int N,Q[maxn];  
22.   
23.   
24. double cross(point a,point b,point c) // c在直线ab右边返回<0  
25. {return (c.y-a.y)*(b.x-a.x)-(b.y-a.y)*(c.x-a.x);}  
26.   
27.   
28. bool operator < (line a,line b) // 直线se右边为可行域  
29. {  
30.     if( fabs(a.k-b.k)<eps )  
31.         return cross(b.s,b.e,a.s)<0;  
32.     return a.k<b.k;  
33. }  
34.   
35.   
36. point intersection(point u1,point u2,point v1,point v2){  
37.     point ret=u1;  
38.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
39.         /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
40.     ret.x+=(u2.x-u1.x)*t;  
41.     ret.y+=(u2.y-u1.y)*t;  
42.     return ret;  
43. }  
44.   
45.   
46. double HalfInSec()  
47. {  
48.     int i,j,k,l;  
49.     sort(L,L+N);  // 极角[-pi,pi]排序  
50.       
51.     for(i=1,j=0;i<N;i++)  // 去掉相同极角直线  
52.         if( fabs(L[i].k-L[j].k)>eps )   
53.         L[++j] = L[i];  
54.     N = j+1;  
55.   
56.   
57.     k = 0,l = 1;  
58.     Q[0] = 0,Q[1] = 1;  
59.     S[1] = intersection(L[0].s,L[0].e,L[1].s,L[1].e);  
60.     for(i=2;i<N;i++)  
61.     {  
62.         while( k<l && cross(L[i].s,L[i].e,S[l])>eps )   
63.             l--;  
64.         while( k<l && cross(L[i].s,L[i].e,S[k+1])>eps )   
65.             k++;  
66.         Q[++l] = i;  
67.         S[l] = intersection(L[Q[l-1]].s,L[Q[l-1]].e,L[i].s,L[i].e);  
68.     }  
69.   
70.   
71.     while( k<l && cross(L[Q[k]].s,L[Q[k]].e,S[l])>eps )  
72.         l--;  
73.     while( k<l && cross(L[Q[l]].s,L[Q[l]].e,S[k+1])>eps )  
74.         k++;  
75.     S[k] = intersection(L[Q[l]].s,L[Q[l]].e,L[Q[k]].s,L[Q[k]].e);  
76.     S[++l] = S[k];  
77.   
78.   
79.     double s = 0;  
80.     for(i=k;i<l;i++)  
81.         s += S[i].y*S[i+1].x-S[i+1].y*S[i].x;  
82.     return fabs(s/2);  
83. }  
84. int main()  
85. {  
86.     int i,j,k,l;  
87.     scanf("%d",&N);  
88.     for(i=0;i<N;i++)  
89.         scanf("%lf%lf%lf%lf",&L[i].e.x,&L[i].e.y,&L[i].s.x,&L[i].s.y);  
90.     L[N].s.x = 0,L[N].s.y = 0;  
91.     L[N+1].s.x = 10000,L[N+1].s.y = 0;  
92.     L[N+2].s.x = 10000,L[N+2].s.y = 10000;  
93.     L[N+3].s.x = 0,L[N+3].s.y = 10000;  
94.     L[N].e = L[N+3].s;  
95.     L[N+1].e = L[N].s;  
96.     L[N+2].e = L[N+1].s;  
97.     L[N+3].e = L[N+2].s;  
98.     N += 4;  
99.     for(i=0;i<N;i++)  
100.         L[i].k = atan2(L[i].s.y-L[i].e.y,L[i].s.x-L[i].e.x);  
101.     printf("%.1lf\n",HalfInSec());  
102. }  

[cpp]  view plain  copy

1. //N^2  
2. #include <iostream>  
3. #include <math.h>  
4. #include <string>  
5. using namespace std;  
6. #define maxn 2005  
7. int N;  
8. struct point  
9. {  
10.     double x, y;  
11. }P[maxn];  
12.   
13.   
14. double xmul(point a,point b,point c)  
15. {return (c.y-a.y)*(b.x-a.x)-(b.y-a.y)*(c.x-a.x);}  
16.   
17.   
18. point intersection(point u1,point u2,point v1,point v2){  
19.     point ret=u1;  
20.     double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))  
21.         /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));  
22.     ret.x+=(u2.x-u1.x)*t;  
23.     ret.y+=(u2.y-u1.y)*t;  
24.     return ret;  
25. }  
26.   
27.   
28. int HalfInteSec(double move)  
29. {  
30.     int i,j;  
31.     point tt[maxn],tp[maxn],ns,ne;  
32.     int len,tlen;  
33.     len = N;  
34.     memcpy(tp,P,sizeof(P));  
35.     for(i=0;i<N;i++)  
36.     {  
37.         point s,e;  
38.         double dx,dy;  
39.           
40.           
41.         dx = P[i].y-P[i+1].y;dy = P[i+1].x-P[i].x;  
42.         dx = P[i].y-P[i+1].y;dy = P[i+1].x-P[i].x;  
43.   
44.   
45.         s.x = P[i+1].x+dx*(move/sqrt(dx*dx+dy*dy));  
46.         s.y = P[i+1].y+dy*(move/sqrt(dx*dx+dy*dy));  
47.           
48.         e.x = P[i].x+dx*(move/sqrt(dx*dx+dy*dy));  
49.         e.y = P[i].y+dy*(move/sqrt(dx*dx+dy*dy));  
50.   
51.   
52.   
53.   
54.         tlen = 0;  
55.         for(j=0;j<len;j++)  
56.         {  
57.             ns = tp[j];    ne = tp[j+1];  
58.               
59.             if( xmul(s,e,ns)<=0 )  
60.                 tt[tlen++] = ns;  
61.               
62.             if( xmul(s,e,ns)*xmul(s,e,ne)<0 )  
63.                 tt[tlen++] = intersection(s,e,ns,ne);  
64.         }  
65.         tt[tlen] = tt[0];  
66.         memcpy(tp,tt,sizeof(tt));  
67.         len = tlen;  
68.     }  
69.     return len;   
70. }  
71.   
72.   
73. int main()  
74. {  
75.     int i,j,k,l;  
76.     while( scanf("%d",&N),N )  
77.     {  
78.         for(i=0;i<N;i++)  
79.             scanf("%lf%lf",&P[i].x,&P[i].y);  
80.         P[N] = P[0];  
81.         double min = 0,max = 10000,mid;  
82.         while( max-min>0.0000001 )  
83.         {  
84.             mid = (min+max)/2;  
85.             l = HalfInteSec(mid);  
86.             if( l==0 )  
87.                 max = mid;  
88.             else  
89.                 min = mid;  
90.         }  
91.         printf("%.6lf\n",min);  
92.     }  
93. }