Yangk's-计算几何
码 计算几何模板
原博客在这
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 多边形
#include <stdlib.h> #include <math.h> #define MAXN 1000 #define offset 10000 #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))<eps) #define _sign(x) ((x)>eps?1:((x)<-eps?2:0)) struct point{double x,y;}; struct line{point a,b;}; double xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } //判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线 int is_convex(int n,point* p){ int i,s[3]={1,1,1}; for (i=0;i<n&&s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0; return s[1]|s[2]; } //判定凸多边形,顶点按顺时针或逆时针给出,不允许相邻边共线 int is_convex_v2(int n,point* p){ int i,s[3]={1,1,1}; for (i=0;i<n&&s[0]&&s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))]=0; return s[0]&&s[1]|s[2]; } //判点在凸多边形内或多边形边上,顶点按顺时针或逆时针给出 int inside_convex(point q,int n,point* p){ int i,s[3]={1,1,1}; for (i=0;i<n&&s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0; return s[1]|s[2]; } //判点在凸多边形内,顶点按顺时针或逆时针给出,在多边形边上返回0 int inside_convex_v2(point q,int n,point* p){ int i,s[3]={1,1,1}; for (i=0;i<n&&s[0]&&s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],q,p[i]))]=0; return s[0]&&s[1]|s[2]; } //判点在任意多边形内,顶点按顺时针或逆时针给出 //on_edge表示点在多边形边上时的返回值,offset为多边形坐标上限 int inside_polygon(point q,int n,point* p,int on_edge=1){ point q2; int i=0,count; while (i<n) for (count=i=0,q2.x=rand()+offset,q2.y=rand()+offset;i<n;i++) 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) return on_edge; else if (zero(xmult(q,q2,p[i]))) break; 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) count++; return count&1; } inline int opposite_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps; } inline int dot_online_in(point p,point l1,point l2){ return zero(xmult(p,l1,l2))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps; } //判线段在任意多边形内,顶点按顺时针或逆时针给出,与边界相交返回1 int inside_polygon(point l1,point l2,int n,point* p){ point t[MAXN],tt; int i,j,k=0; if (!inside_polygon(l1,n,p)||!inside_polygon(l2,n,p)) return 0; for (i=0;i<n;i++) if (opposite_side(l1,l2,p[i],p[(i+1)%n])&&opposite_side(p[i],p[(i+1)%n],l1,l2)) return 0; else if (dot_online_in(l1,p[i],p[(i+1)%n])) t[k++]=l1; else if (dot_online_in(l2,p[i],p[(i+1)%n])) t[k++]=l2; else if (dot_online_in(p[i],l1,l2)) t[k++]=p[i]; for (i=0;i<k;i++) for (j=i+1;j<k;j++){ tt.x=(t[i].x+t[j].x)/2; tt.y=(t[i].y+t[j].y)/2; if (!inside_polygon(tt,n,p)) return 0; } return 1; } point intersection(line u,line v){ point ret=u.a; 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)) /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x)); ret.x+=(u.b.x-u.a.x)*t; ret.y+=(u.b.y-u.a.y)*t; return ret; } point barycenter(point a,point b,point c){ line u,v; u.a.x=(a.x+b.x)/2; u.a.y=(a.y+b.y)/2; u.b=c; v.a.x=(a.x+c.x)/2; v.a.y=(a.y+c.y)/2; v.b=b; return intersection(u,v); } //多边形重心 point barycenter(int n,point* p){ point ret,t; double t1=0,t2; int i; ret.x=ret.y=0; for (i=1;i<n-1;i++) if (fabs(t2=xmult(p[0],p[i],p[i+1]))>eps){ t=barycenter(p[0],p[i],p[i+1]); ret.x+=t.x*t2; ret.y+=t.y*t2; t1+=t2; } if (fabs(t1)>eps) ret.x/=t1,ret.y/=t1; return ret; }
1.4多边形切割
//多边形切割 //可用于半平面交 #define MAXN 100 #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))<eps) struct point{double x,y;}; double xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } int same_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps; } point intersection(point u1,point u2,point v1,point v2){ point ret=u1; double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)) /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x)); ret.x+=(u2.x-u1.x)*t; ret.y+=(u2.y-u1.y)*t; return ret; } //将多边形沿l1,l2确定的直线切割在side侧切割,保证l1,l2,side不共线 void polygon_cut(int& n,point* p,point l1,point l2,point side){ point pp[MAXN]; int m=0,i; for (i=0;i<n;i++){ if (same_side(p[i],side,l1,l2)) pp[m++]=p[i]; 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)))) pp[m++]=intersection(p[i],p[(i+1)%n],l1,l2); } for (n=i=0;i<m;i++) if (!i||!zero(pp[i].x-pp[i-1].x)||!zero(pp[i].y-pp[i-1].y)) p[n++]=pp[i]; if (zero(p[n-1].x-p[0].x)&&zero(p[n-1].y-p[0].y)) n--; if (n<3) n=0; }
1.5 浮点函数
//浮点几何函数库 #include <math.h> #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))<eps) struct point{double x,y;}; struct line{point a,b;}; //计算cross product (P1-P0)x(P2-P0) double xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } double xmult(double x1,double y1,double x2,double y2,double x0,double y0){ return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0); } //计算dot product (P1-P0).(P2-P0) double dmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y); } double dmult(double x1,double y1,double x2,double y2,double x0,double y0){ return (x1-x0)*(x2-x0)+(y1-y0)*(y2-y0); } //两点距离 double distance(point p1,point p2){ return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } double distance(double x1,double y1,double x2,double y2){ return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); } //判三点共线 int dots_inline(point p1,point p2,point p3){ return zero(xmult(p1,p2,p3)); } int dots_inline(double x1,double y1,double x2,double y2,double x3,double y3){ return zero(xmult(x1,y1,x2,y2,x3,y3)); } //判点是否在线段上,包括端点 int dot_online_in(point p,line l){ 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; } int dot_online_in(point p,point l1,point l2){ return zero(xmult(p,l1,l2))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps; } int dot_online_in(double x,double y,double x1,double y1,double x2,double y2){ return zero(xmult(x,y,x1,y1,x2,y2))&&(x1-x)*(x2-x)<eps&&(y1-y)*(y2-y)<eps; } //判点是否在线段上,不包括端点 int dot_online_ex(point p,line l){ 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)); } int dot_online_ex(point p,point l1,point l2){ 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)); } int dot_online_ex(double x,double y,double x1,double y1,double x2,double y2){ return dot_online_in(x,y,x1,y1,x2,y2)&&(!zero(x-x1)||!zero(y-y1))&&(!zero(x-x2)||!zero(y-y2)); } //判两点在线段同侧,点在线段上返回0 int same_side(point p1,point p2,line l){ return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)>eps; } int same_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps; } //判两点在线段异侧,点在线段上返回0 int opposite_side(point p1,point p2,line l){ return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps; } int opposite_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps; } //判两直线平行 int parallel(line u,line v){ 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)); } int parallel(point u1,point u2,point v1,point v2){ return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y)); } //判两直线垂直 int perpendicular(line u,line v){ 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)); } int perpendicular(point u1,point u2,point v1,point v2){ return zero((u1.x-u2.x)*(v1.x-v2.x)+(u1.y-u2.y)*(v1.y-v2.y)); } //判两线段相交,包括端点和部分重合 int intersect_in(line u,line v){ if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b)) return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u); 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); } int intersect_in(point u1,point u2,point v1,point v2){ if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2)) return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2); 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); } //判两线段相交,不包括端点和部分重合 int intersect_ex(line u,line v){ return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u); } int intersect_ex(point u1,point u2,point v1,point v2){ return opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2); } //计算两直线交点,注意事先判断直线是否平行! //线段交点请另外判线段相交(同时还是要判断是否平行!) point intersection(line u,line v){ point ret=u.a; 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)) /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x)); ret.x+=(u.b.x-u.a.x)*t; ret.y+=(u.b.y-u.a.y)*t; return ret; } point intersection(point u1,point u2,point v1,point v2){ point ret=u1; double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)) /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x)); ret.x+=(u2.x-u1.x)*t; ret.y+=(u2.y-u1.y)*t; return ret; } //点到直线上的最近点 point ptoline(point p,line l){ point t=p; t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x; return intersection(p,t,l.a,l.b); } point ptoline(point p,point l1,point l2){ point t=p; t.x+=l1.y-l2.y,t.y+=l2.x-l1.x; return intersection(p,t,l1,l2); } //点到直线距离 double disptoline(point p,line l){ return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b); } double disptoline(point p,point l1,point l2){ return fabs(xmult(p,l1,l2))/distance(l1,l2); } double disptoline(double x,double y,double x1,double y1,double x2,double y2){ return fabs(xmult(x,y,x1,y1,x2,y2))/distance(x1,y1,x2,y2); } //点到线段上的最近点 point ptoseg(point p,line l){ point t=p; t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x; if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps) return distance(p,l.a)<distance(p,l.b)?l.a:l.b; return intersection(p,t,l.a,l.b); } point ptoseg(point p,point l1,point l2){ point t=p; t.x+=l1.y-l2.y,t.y+=l2.x-l1.x; if (xmult(l1,t,p)*xmult(l2,t,p)>eps) return distance(p,l1)<distance(p,l2)?l1:l2; return intersection(p,t,l1,l2); } //点到线段距离 double disptoseg(point p,line l){ point t=p; t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x; if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps) return distance(p,l.a)<distance(p,l.b)?distance(p,l.a):distance(p,l.b); return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b); } double disptoseg(point p,point l1,point l2){ point t=p; t.x+=l1.y-l2.y,t.y+=l2.x-l1.x; if (xmult(l1,t,p)*xmult(l2,t,p)>eps) return distance(p,l1)<distance(p,l2)?distance(p,l1):distance(p,l2); return fabs(xmult(p,l1,l2))/distance(l1,l2); } //矢量V以P为顶点逆时针旋转angle并放大scale倍 point rotate(point v,point p,double angle,double scale){ point ret=p; v.x-=p.x,v.y-=p.y; p.x=scale*cos(angle); p.y=scale*sin(angle); ret.x+=v.x*p.x-v.y*p.y; ret.y+=v.x*p.y+v.y*p.x; return ret; } //p点关于直线L的对称点 ponit symmetricalPointofLine(point p, line L) { point p2; double d; d = L.a * L.a + L.b * L.b; p2.x = (L.b * L.b * p.x - L.a * L.a * p.x - 2 * L.a * L.b * p.y - 2 * L.a * L.c) / d; p2.y = (L.a * L.a * p.y - L.b * L.b * p.y - 2 * L.a * L.b * p.x - 2 * L.b * L.c) / d; return p2; } //求两点的平分线 line bisector(point& a, point& b) { line ab, ans; ab.set(a, b); double midx = (a.x + b.x)/2.0, midy = (a.y + b.y)/2.0; ans.a = -ab.b, ans.b = -ab.a, ans.c = -ab.b * midx + ab.a * midy; return ans; } // 已知入射线、镜面，求反射线。 // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数; a2,b2,c2为入射光直线方程系数; a,b,c为反射光直线方程系数. // 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>. // 不要忘记结果中可能会有"negative zeros" void reflect(double a1,double b1,double c1, double a2,double b2,double c2, double &a,double &b,double &c) { double n,m; double tpb,tpa; tpb=b1*b2+a1*a2; tpa=a2*b1-a1*b2; m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1); n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1); if(fabs(a1*b2-a2*b1)<1e-20) { a=a2;b=b2;c=c2; return; } double xx,yy; //(xx,yy)是入射线与镜面的交点。 xx=(b1*c2-b2*c1)/(a1*b2-a2*b1); yy=(a2*c1-a1*c2)/(a1*b2-a2*b1); a=n; b=-m; c=m*yy-xx*n; }
1.6 面积
#include <math.h> struct point{double x,y;}; //计算cross product (P1-P0)x(P2-P0) double xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } double xmult(double x1,double y1,double x2,double y2,double x0,double y0){ return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0); } //计算三角形面积,输入三顶点 double area_triangle(point p1,point p2,point p3){ return fabs(xmult(p1,p2,p3))/2; } double area_triangle(double x1,double y1,double x2,double y2,double x3,double y3){ return fabs(xmult(x1,y1,x2,y2,x3,y3))/2; } //计算三角形面积,输入三边长 double area_triangle(double a,double b,double c){ double s=(a+b+c)/2; return sqrt(s*(s-a)*(s-b)*(s-c)); } //计算多边形面积,顶点按顺时针或逆时针给出 double area_polygon(int n,point* p){ double s1=0,s2=0; int i; for (i=0;i<n;i++) s1+=p[(i+1)%n].y*p[i].x,s2+=p[(i+1)%n].y*p[(i+2)%n].x; return fabs(s1-s2)/2; }
1.7球面
#include <math.h> const double pi=acos(-1); //计算圆心角lat表示纬度,-90<=w<=90,lng表示经度 //返回两点所在大圆劣弧对应圆心角,0<=angle<=pi double angle(double lng1,double lat1,double lng2,double lat2){ double dlng=fabs(lng1-lng2)*pi/180; while (dlng>=pi+pi) dlng-=pi+pi; if (dlng>pi) dlng=pi+pi-dlng; lat1*=pi/180,lat2*=pi/180; return acos(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2)); } //计算距离,r为球半径 double line_dist(double r,double lng1,double lat1,double lng2,double lat2){ double dlng=fabs(lng1-lng2)*pi/180; while (dlng>=pi+pi) dlng-=pi+pi; if (dlng>pi) dlng=pi+pi-dlng; lat1*=pi/180,lat2*=pi/180; return r*sqrt(2-2*(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2))); } //计算球面距离,r为球半径 inline double sphere_dist(double r,double lng1,double lat1,double lng2,double lat2){ return r*angle(lng1,lat1,lng2,lat2); }
1.8三角形
#include <math.h> struct point{double x,y;}; struct line{point a,b;}; double distance(point p1,point p2){ return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } point intersection(line u,line v){ point ret=u.a; 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)) /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x)); ret.x+=(u.b.x-u.a.x)*t; ret.y+=(u.b.y-u.a.y)*t; return ret; } //外心 point circumcenter(point a,point b,point c){ line u,v; u.a.x=(a.x+b.x)/2; u.a.y=(a.y+b.y)/2; u.b.x=u.a.x-a.y+b.y; u.b.y=u.a.y+a.x-b.x; v.a.x=(a.x+c.x)/2; v.a.y=(a.y+c.y)/2; v.b.x=v.a.x-a.y+c.y; v.b.y=v.a.y+a.x-c.x; return intersection(u,v); } //内心 point incenter(point a,point b,point c){ line u,v; double m,n; u.a=a; m=atan2(b.y-a.y,b.x-a.x); n=atan2(c.y-a.y,c.x-a.x); u.b.x=u.a.x+cos((m+n)/2); u.b.y=u.a.y+sin((m+n)/2); v.a=b; m=atan2(a.y-b.y,a.x-b.x); n=atan2(c.y-b.y,c.x-b.x); v.b.x=v.a.x+cos((m+n)/2); v.b.y=v.a.y+sin((m+n)/2); return intersection(u,v); } //垂心 point perpencenter(point a,point b,point c){ line u,v; u.a=c; u.b.x=u.a.x-a.y+b.y; u.b.y=u.a.y+a.x-b.x; v.a=b; v.b.x=v.a.x-a.y+c.y; v.b.y=v.a.y+a.x-c.x; return intersection(u,v); } //重心 //到三角形三顶点距离的平方和最小的点 //三角形内到三边距离之积最大的点 point barycenter(point a,point b,point c){ line u,v; u.a.x=(a.x+b.x)/2; u.a.y=(a.y+b.y)/2; u.b=c; v.a.x=(a.x+c.x)/2; v.a.y=(a.y+c.y)/2; v.b=b; return intersection(u,v); } //费马点 //到三角形三顶点距离之和最小的点 point fermentpoint(point a,point b,point c){ point u,v; double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y); int i,j,k; u.x=(a.x+b.x+c.x)/3; u.y=(a.y+b.y+c.y)/3; while (step>1e-10) for (k=0;k<10;step/=2,k++) for (i=-1;i<=1;i++) for (j=-1;j<=1;j++){ v.x=u.x+step*i; v.y=u.y+step*j; if (distance(u,a)+distance(u,b)+distance(u,c)>distance(v,a)+distance(v,b)+distance(v,c)) u=v; } return u; } //求曲率半径 三角形内最大可围成面积 #include<iostream> #include<cmath> using namespace std; const double pi=3.14159265358979; int main() { double a,b,c,d,p,s,r,ans,R,x,l; int T=0; while(cin>>a>>b>>c>>d&&a+b+c+d) { T++; l=a+b+c; p=l/2; s=sqrt(p*(p-a)*(p-b)*(p-c)); R= s /p; if (d >= l) ans = s; else if(2*pi*R>=d) ans=d*d/(4*pi); else { r = (l-d)/((l/R)-(2*pi)); x = r*r*s/(R*R); ans = s - x + pi * r * r; } printf("Case %d: %.2lf\n",T,ans); } return 0; }
1.9三维几何
//三维几何函数库 #include <math.h> #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))<eps) struct point3{double x,y,z;}; struct line3{point3 a,b;}; struct plane3{point3 a,b,c;}; //计算cross product U x V point3 xmult(point3 u,point3 v){ point3 ret; ret.x=u.y*v.z-v.y*u.z; ret.y=u.z*v.x-u.x*v.z; ret.z=u.x*v.y-u.y*v.x; return ret; } //计算dot product U . V double dmult(point3 u,point3 v){ return u.x*v.x+u.y*v.y+u.z*v.z; } //矢量差 U - V point3 subt(point3 u,point3 v){ point3 ret; ret.x=u.x-v.x; ret.y=u.y-v.y; ret.z=u.z-v.z; return ret; } //取平面法向量 point3 pvec(plane3 s){ return xmult(subt(s.a,s.b),subt(s.b,s.c)); } point3 pvec(point3 s1,point3 s2,point3 s3){ return xmult(subt(s1,s2),subt(s2,s3)); } //两点距离,单参数取向量大小 double distance(point3 p1,point3 p2){ 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)); } //向量大小 double vlen(point3 p){ return sqrt(p.x*p.x+p.y*p.y+p.z*p.z); } //判三点共线 int dots_inline(point3 p1,point3 p2,point3 p3){ return vlen(xmult(subt(p1,p2),subt(p2,p3)))<eps; } //判四点共面 int dots_onplane(point3 a,point3 b,point3 c,point3 d){ return zero(dmult(pvec(a,b,c),subt(d,a))); } //判点是否在线段上,包括端点和共线 int dot_online_in(point3 p,line3 l){ return zero(vlen(xmult(subt(p,l.a),subt(p,l.b))))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&& (l.a.y-p.y)*(l.b.y-p.y)<eps&&(l.a.z-p.z)*(l.b.z-p.z)<eps; } int dot_online_in(point3 p,point3 l1,point3 l2){ return zero(vlen(xmult(subt(p,l1),subt(p,l2))))&&(l1.x-p.x)*(l2.x-p.x)<eps&& (l1.y-p.y)*(l2.y-p.y)<eps&&(l1.z-p.z)*(l2.z-p.z)<eps; } //判点是否在线段上,不包括端点 int dot_online_ex(point3 p,line3 l){ return dot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y)||!zero(p.z-l.a.z))&& (!zero(p.x-l.b.x)||!zero(p.y-l.b.y)||!zero(p.z-l.b.z)); } int dot_online_ex(point3 p,point3 l1,point3 l2){ return dot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y)||!zero(p.z-l1.z))&& (!zero(p.x-l2.x)||!zero(p.y-l2.y)||!zero(p.z-l2.z)); } //判点是否在空间三角形上,包括边界,三点共线无意义 int dot_inplane_in(point3 p,plane3 s){ return zero(vlen(xmult(subt(s.a,s.b),subt(s.a,s.c)))-vlen(xmult(subt(p,s.a),subt(p,s.b)))- vlen(xmult(subt(p,s.b),subt(p,s.c)))-vlen(xmult(subt(p,s.c),subt(p,s.a)))); } int dot_inplane_in(point3 p,point3 s1,point3 s2,point3 s3){ return zero(vlen(xmult(subt(s1,s2),subt(s1,s3)))-vlen(xmult(subt(p,s1),subt(p,s2)))- vlen(xmult(subt(p,s2),subt(p,s3)))-vlen(xmult(subt(p,s3),subt(p,s1)))); } //判点是否在空间三角形上,不包括边界,三点共线无意义 int dot_inplane_ex(point3 p,plane3 s){ return dot_inplane_in(p,s)&&vlen(xmult(subt(p,s.a),subt(p,s.b)))>eps&& vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps; } int dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3){ return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&& vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps; } //判两点在线段同侧,点在线段上返回0,不共面无意义 int same_side(point3 p1,point3 p2,line3 l){ return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))>eps; } int same_side(point3 p1,point3 p2,point3 l1,point3 l2){ return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))>eps; } //判两点在线段异侧,点在线段上返回0,不共面无意义 int opposite_side(point3 p1,point3 p2,line3 l){ return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))<-eps; } int opposite_side(point3 p1,point3 p2,point3 l1,point3 l2){ return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))<-eps; } //判两点在平面同侧,点在平面上返回0 int same_side(point3 p1,point3 p2,plane3 s){ return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps; } int same_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){ return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))>eps; } //判两点在平面异侧,点在平面上返回0 int opposite_side(point3 p1,point3 p2,plane3 s){ return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps; } int opposite_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){ return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))<-eps; } //判两直线平行 int parallel(line3 u,line3 v){ return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))<eps; } int parallel(point3 u1,point3 u2,point3 v1,point3 v2){ return vlen(xmult(subt(u1,u2),subt(v1,v2)))<eps; } //判两平面平行 int parallel(plane3 u,plane3 v){ return vlen(xmult(pvec(u),pvec(v)))<eps; } int parallel(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){ return vlen(xmult(pvec(u1,u2,u3),pvec(v1,v2,v3)))<eps; } //判直线与平面平行 int parallel(line3 l,plane3 s){ return zero(dmult(subt(l.a,l.b),pvec(s))); } int parallel(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ return zero(dmult(subt(l1,l2),pvec(s1,s2,s3))); } //判两直线垂直 int perpendicular(line3 u,line3 v){ return zero(dmult(subt(u.a,u.b),subt(v.a,v.b))); } int perpendicular(point3 u1,point3 u2,point3 v1,point3 v2){ return zero(dmult(subt(u1,u2),subt(v1,v2))); } //判两平面垂直 int perpendicular(plane3 u,plane3 v){ return zero(dmult(pvec(u),pvec(v))); } int perpendicular(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){ return zero(dmult(pvec(u1,u2,u3),pvec(v1,v2,v3))); } //判直线与平面平行 int perpendicular(line3 l,plane3 s){ return vlen(xmult(subt(l.a,l.b),pvec(s)))<eps; } int perpendicular(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ return vlen(xmult(subt(l1,l2),pvec(s1,s2,s3)))<eps; } //判两线段相交,包括端点和部分重合 int intersect_in(line3 u,line3 v){ if (!dots_onplane(u.a,u.b,v.a,v.b)) return 0; if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b)) return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u); 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); } int intersect_in(point3 u1,point3 u2,point3 v1,point3 v2){ if (!dots_onplane(u1,u2,v1,v2)) return 0; if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2)) return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2); 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); } //判两线段相交,不包括端点和部分重合 int intersect_ex(line3 u,line3 v){ 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); } int intersect_ex(point3 u1,point3 u2,point3 v1,point3 v2){ return dots_onplane(u1,u2,v1,v2)&&opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2); } //判线段与空间三角形相交,包括交于边界和(部分)包含 int intersect_in(line3 l,plane3 s){ return !same_side(l.a,l.b,s)&&!same_side(s.a,s.b,l.a,l.b,s.c)&& !same_side(s.b,s.c,l.a,l.b,s.a)&&!same_side(s.c,s.a,l.a,l.b,s.b); } int intersect_in(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ return !same_side(l1,l2,s1,s2,s3)&&!same_side(s1,s2,l1,l2,s3)&& !same_side(s2,s3,l1,l2,s1)&&!same_side(s3,s1,l1,l2,s2); } //判线段与空间三角形相交,不包括交于边界和(部分)包含 int intersect_ex(line3 l,plane3 s){ return opposite_side(l.a,l.b,s)&&opposite_side(s.a,s.b,l.a,l.b,s.c)&& opposite_side(s.b,s.c,l.a,l.b,s.a)&&opposite_side(s.c,s.a,l.a,l.b,s.b); } int intersect_ex(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ return opposite_side(l1,l2,s1,s2,s3)&&opposite_side(s1,s2,l1,l2,s3)&& opposite_side(s2,s3,l1,l2,s1)&&opposite_side(s3,s1,l1,l2,s2); } //计算两直线交点,注意事先判断直线是否共面和平行! //线段交点请另外判线段相交(同时还是要判断是否平行!) point3 intersection(line3 u,line3 v){ point3 ret=u.a; 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)) /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x)); ret.x+=(u.b.x-u.a.x)*t; ret.y+=(u.b.y-u.a.y)*t; ret.z+=(u.b.z-u.a.z)*t; return ret; } point3 intersection(point3 u1,point3 u2,point3 v1,point3 v2){ point3 ret=u1; double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)) /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x)); ret.x+=(u2.x-u1.x)*t; ret.y+=(u2.y-u1.y)*t; ret.z+=(u2.z-u1.z)*t; return ret; } //计算直线与平面交点,注意事先判断是否平行,并保证三点不共线! //线段和空间三角形交点请另外判断 point3 intersection(line3 l,plane3 s){ point3 ret=pvec(s); 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))/ (ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z)); ret.x=l.a.x+(l.b.x-l.a.x)*t; ret.y=l.a.y+(l.b.y-l.a.y)*t; ret.z=l.a.z+(l.b.z-l.a.z)*t; return ret; } point3 intersection(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ point3 ret=pvec(s1,s2,s3); double t=(ret.x*(s1.x-l1.x)+ret.y*(s1.y-l1.y)+ret.z*(s1.z-l1.z))/ (ret.x*(l2.x-l1.x)+ret.y*(l2.y-l1.y)+ret.z*(l2.z-l1.z)); ret.x=l1.x+(l2.x-l1.x)*t; ret.y=l1.y+(l2.y-l1.y)*t; ret.z=l1.z+(l2.z-l1.z)*t; return ret; } //计算两平面交线,注意事先判断是否平行,并保证三点不共线! line3 intersection(plane3 u,plane3 v){ line3 ret; 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); 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); return ret; } line3 intersection(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){ line3 ret; ret.a=parallel(v1,v2,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v1,v2,u1,u2,u3); ret.b=parallel(v3,v1,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v3,v1,u1,u2,u3); return ret; } //点到直线距离 double ptoline(point3 p,line3 l){ return vlen(xmult(subt(p,l.a),subt(l.b,l.a)))/distance(l.a,l.b); } double ptoline(point3 p,point3 l1,point3 l2){ return vlen(xmult(subt(p,l1),subt(l2,l1)))/distance(l1,l2); } //点到平面距离 double ptoplane(point3 p,plane3 s){ return fabs(dmult(pvec(s),subt(p,s.a)))/vlen(pvec(s)); } double ptoplane(point3 p,point3 s1,point3 s2,point3 s3){ return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3)); } //直线到直线距离 double linetoline(line3 u,line3 v){ point3 n=xmult(subt(u.a,u.b),subt(v.a,v.b)); return fabs(dmult(subt(u.a,v.a),n))/vlen(n); } double linetoline(point3 u1,point3 u2,point3 v1,point3 v2){ point3 n=xmult(subt(u1,u2),subt(v1,v2)); return fabs(dmult(subt(u1,v1),n))/vlen(n); } //两直线夹角cos值 double angle_cos(line3 u,line3 v){ return dmult(subt(u.a,u.b),subt(v.a,v.b))/vlen(subt(u.a,u.b))/vlen(subt(v.a,v.b)); } double angle_cos(point3 u1,point3 u2,point3 v1,point3 v2){ return dmult(subt(u1,u2),subt(v1,v2))/vlen(subt(u1,u2))/vlen(subt(v1,v2)); } //两平面夹角cos值 double angle_cos(plane3 u,plane3 v){ return dmult(pvec(u),pvec(v))/vlen(pvec(u))/vlen(pvec(v)); } double angle_cos(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){ return dmult(pvec(u1,u2,u3),pvec(v1,v2,v3))/vlen(pvec(u1,u2,u3))/vlen(pvec(v1,v2,v3)); } //直线平面夹角sin值 double angle_sin(line3 l,plane3 s){ return dmult(subt(l.a,l.b),pvec(s))/vlen(subt(l.a,l.b))/vlen(pvec(s)); } double angle_sin(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){ return dmult(subt(l1,l2),pvec(s1,s2,s3))/vlen(subt(l1,l2))/vlen(pvec(s1,s2,s3)); }
1.10 凸包
//水平序 #define maxn 100005 struct point {double x,y;}p[maxn],s[maxn]; bool operator < (point a,point b) {return a.x<b.x || a.x==b.x&&a.y<b.y;} int n,f; double cp(point a,point b,point c) {return (c.y-a.y)*(b.x-a.x)-(b.y-a.y)*(c.x-a.x);} void Convex(point *p,int &n) { sort(p,p+n); int i,j,r,top,m; s[0] = p[0];s[1] = p[1];top = 1; for(i=2;i<n;i++) { while( top>0 && cp(p[i],s[top],s[top-1])>=0 ) top--; top++;s[top] = p[i]; } m = top; top++;s[top] = p[n-2]; for(i=n-3;i>=0;i--) { while( top>m && cp(p[i],s[top],s[top-1])>=0 ) top--; top++;s[top] = p[i]; } top--; n = top+1; } 极角序 #include <stdio.h> #include <string.h> #include <algorithm> #include <math.h> using namespace std; #define maxn 100005 int N; struct A { int x,y; int v,l; }P[maxn]; int xmult(int x1,int y1,int x2,int y2,int x3,int y3) { return (y2-y1)*(x3-x1)-(y3-y1)*(x2-x1); } void swap(A &a,A &b) { A t = a;a = b,b = t; } bool operator < (A a,A b) { int k = xmult(P[0].x,P[0].y,a.x,a.y,b.x,b.y); if( k<0 ) return 1; else if( k==0 ) { if( abs(P[0].x-a.x)<abs(P[0].x-b.x) ) return 1; if( abs(P[0].y-a.y)<abs(P[0].y-b.y) ) return 1; } return 0; } void Grem_scan(int n) { int i,j,k,l; k = 0x7fffffff; for(i=0;i<n;i++) if( P[i].x<k || P[i].x==k && P[i].y<P[l].y ) k = P[i].x,l = i; swap(P[l],P[0]); sort(P+1,P+n); l = 3; for(i=3;i<n;i++) { 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 ) l--; P[l++] = P[i]; } } main() { int i,j,k,l; N = 0; while( scanf("%d%d",&P[N].x,&P[N].y)!=EOF ) N++; Grem_scan(N); for(i=0;i<N;i++) if( P[i].x==0 && P[i].y==0 ) break; k = i++; printf("(0,0)\n"); while( i!=k ) printf("(%d,%d)\n",P[i].x,P[i].y),i = (i+1)%N; } //卷包裹法 #include <stdio.h> #include <string.h> #include <algorithm> using namespace std; #define maxn 55 struct A { int x,y; }P[maxn]; int T,N; bool B[maxn]; int as[maxn],L; int xmult(A a,A b,A c) { return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x); } int main() { int i,j,k,l; scanf("%d",&T); while( T-- ) { scanf("%d",&N); k = 0x7ffffff; for(i=0;i<N;i++) { scanf("%d%d%d",&j,&P[i].x,&P[i].y); if( P[i].y<k ) k = P[i].y,l = i; } memset(B,0,sizeof(B)); B[l] = 1; as[0] = l; L = 1; while( 1 ) { A a,b; if( L==1 ) a.x = 0,a.y = P[as[0]].y; else a = P[as[L-2]]; b = P[as[L-1]]; k = -1; for(i=0;i<N;i++) { if( B[i] ) continue; if( xmult(a,b,P[i])<0 ) continue; 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 ) k = i; } if( k==-1 ) break; B[k] = 1; as[L++] = k; } printf("%d ",L); for(i=0;i<L;i++) printf("%d ",as[i]+1); printf("\n"); } }
1.11 网格
#define abs(x) ((x)>0?(x):-(x)) struct point{int x,y;}; int gcd(int a,int b){return b?gcd(b,a%b):a;} //多边形上的网格点个数 int grid_onedge(int n,point* p){ int i,ret=0; for (i=0;i<n;i++) ret+=gcd(abs(p[i].x-p[(i+1)%n].x),abs(p[i].y-p[(i+1)%n].y)); return ret; } //多边形内的网格点个数 int grid_inside(int n,point* p){ int i,ret=0; for (i=0;i<n;i++) ret+=p[(i+1)%n].y*(p[i].x-p[(i+2)%n].x); return (abs(ret)-grid_onedge(n,p))/2+1; }
1.12 圆
#include <math.h> #define eps 1e-8 struct point{double x,y;}; double xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } double distance(point p1,point p2){ return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } double disptoline(point p,point l1,point l2){ return fabs(xmult(p,l1,l2))/distance(l1,l2); } point intersection(point u1,point u2,point v1,point v2){ point ret=u1; double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)) /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x)); ret.x+=(u2.x-u1.x)*t; ret.y+=(u2.y-u1.y)*t; return ret; } //判直线和圆相交,包括相切 int intersect_line_circle(point c,double r,point l1,point l2){ return disptoline(c,l1,l2)<r+eps; } //判线段和圆相交,包括端点和相切 int intersect_seg_circle(point c,double r,point l1,point l2){ double t1=distance(c,l1)-r,t2=distance(c,l2)-r; point t=c; if (t1<eps||t2<eps) return t1>-eps||t2>-eps; t.x+=l1.y-l2.y; t.y+=l2.x-l1.x; return xmult(l1,c,t)*xmult(l2,c,t)<eps&&disptoline(c,l1,l2)-r<eps; } //判圆和圆相交,包括相切 int intersect_circle_circle(point c1,double r1,point c2,double r2){ return distance(c1,c2)<r1+r2+eps&&distance(c1,c2)>fabs(r1-r2)-eps; } //计算圆上到点p最近点,如p与圆心重合,返回p本身 point dot_to_circle(point c,double r,point p){ point u,v; if (distance(p,c)<eps) return p; u.x=c.x+r*fabs(c.x-p.x)/distance(c,p); u.y=c.y+r*fabs(c.y-p.y)/distance(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1); v.x=c.x-r*fabs(c.x-p.x)/distance(c,p); v.y=c.y-r*fabs(c.y-p.y)/distance(c,p)*((c.x-p.x)*(c.y-p.y)<0?-1:1); return distance(u,p)<distance(v,p)?u:v; } //计算直线与圆的交点,保证直线与圆有交点 //计算线段与圆的交点可用这个函数后判点是否在线段上 void intersection_line_circle(point c,double r,point l1,point l2,point& p1,point& p2){ point p=c; double t; p.x+=l1.y-l2.y; p.y+=l2.x-l1.x; p=intersection(p,c,l1,l2); t=sqrt(r*r-distance(p,c)*distance(p,c))/distance(l1,l2); p1.x=p.x+(l2.x-l1.x)*t; p1.y=p.y+(l2.y-l1.y)*t; p2.x=p.x-(l2.x-l1.x)*t; p2.y=p.y-(l2.y-l1.y)*t; } //计算圆与圆的交点,保证圆与圆有交点,圆心不重合 void intersection_circle_circle(point c1,double r1,point c2,double r2,point& p1,point& p2){ point u,v; double t; t=(1+(r1*r1-r2*r2)/distance(c1,c2)/distance(c1,c2))/2; u.x=c1.x+(c2.x-c1.x)*t; u.y=c1.y+(c2.y-c1.y)*t; v.x=u.x+c1.y-c2.y; v.y=u.y-c1.x+c2.x; intersection_line_circle(c1,r1,u,v,p1,p2); } //将向量p逆时针旋转angle角度 Point Rotate(Point p,double angle) { Point res; res.x=p.x*cos(angle)-p.y*sin(angle); res.y=p.x*sin(angle)+p.y*cos(angle); return res; } //求圆外一点对圆(o,r)的两个切点result1和result2 void TangentPoint_PC(Point poi,Point o,double r,Point &result1,Point &result2) { double line=sqrt((poi.x-o.x)*(poi.x-o.x)+(poi.y-o.y)*(poi.y-o.y)); double angle=acos(r/line); Point unitvector,lin; lin.x=poi.x-o.x; lin.y=poi.y-o.y; unitvector.x=lin.x/sqrt(lin.x*lin.x+lin.y*lin.y)*r; unitvector.y=lin.y/sqrt(lin.x*lin.x+lin.y*lin.y)*r; result1=Rotate(unitvector,-angle); result2=Rotate(unitvector,angle); result1.x+=o.x; result1.y+=o.y; result2.x+=o.x; result2.y+=o.y; return; }
1.13 矢量运算求几何模板
#include <iostream> #include <cmath> #include <vector> #include <algorithm> #define MAX_N 100 using namespace std; /////////////////////////////////////////////////////////////////// //常量区 const double INF = 1e10; // 无穷大 const double EPS = 1e-15; // 计算精度 const int LEFT = 0; // 点在直线左边 const int RIGHT = 1; // 点在直线右边 const int ONLINE = 2; // 点在直线上 const int CROSS = 0; // 两直线相交 const int COLINE = 1; // 两直线共线 const int PARALLEL = 2; // 两直线平行 const int NOTCOPLANAR = 3; // 两直线不共面 const int INSIDE = 1; // 点在图形内部 const int OUTSIDE = 2; // 点在图形外部 const int BORDER = 3; // 点在图形边界 const int BAOHAN = 1; // 大圆包含小圆 const int NEIQIE = 2; // 内切 const int XIANJIAO = 3; // 相交 const int WAIQIE = 4; // 外切 const int XIANLI = 5; // 相离 const double pi = acos(-1.0) //圆周率 /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// //类型定义区 struct Point { // 二维点或矢量 double x, y; double angle, dis; Point() {} Point(double x0, double y0): x(x0), y(y0) {} }; struct Point3D { //三维点或矢量 double x, y, z; Point3D() {} Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {} }; struct Line { // 二维的直线或线段 Point p1, p2; Line() {} Line(Point p10, Point p20): p1(p10), p2(p20) {} }; struct Line3D { // 三维的直线或线段 Point3D p1, p2; Line3D() {} Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {} }; struct Rect { // 用长宽表示矩形的方法 w, h分别表示宽度和高度 double w, h; Rect() {} Rect(double _w,double _h) : w(_w),h(_h) {} }; struct Rect_2 { // 表示矩形，左下角坐标是(xl, yl)，右上角坐标是(xh, yh) double xl, yl, xh, yh; Rect_2() {} Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {} }; struct Circle { //圆 Point c; double r; Circle() {} Circle(Point _c,double _r) :c(_c),r(_r) {} }; typedef vector<Point> Polygon; // 二维多边形 typedef vector<Point> Points; // 二维点集 typedef vector<Point3D> Points3D; // 三维点集 /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// //基本函数区 inline double max(double x,double y) { return x > y ? x : y; } inline double min(double x, double y) { return x > y ? y : x; } inline bool ZERO(double x) // x == 0 { return (fabs(x) < EPS); } inline bool ZERO(Point p) // p == 0 { return (ZERO(p.x) && ZERO(p.y)); } inline bool ZERO(Point3D p) // p == 0 { return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z)); } inline bool EQ(double x, double y) // eqaul, x == y { return (fabs(x - y) < EPS); } inline bool NEQ(double x, double y) // not equal, x != y { return (fabs(x - y) >= EPS); } inline bool LT(double x, double y) // less than, x < y { return ( NEQ(x, y) && (x < y) ); } inline bool GT(double x, double y) // greater than, x > y { return ( NEQ(x, y) && (x > y) ); } inline bool LEQ(double x, double y) // less equal, x <= y { return ( EQ(x, y) || (x < y) ); } inline bool GEQ(double x, double y) // greater equal, x >= y { return ( EQ(x, y) || (x > y) ); } // 注意！！！ // 如果是一个很小的负的浮点数 // 保留有效位数输出的时候会出现-0.000这样的形式， // 前面多了一个负号 // 这就会导致错误！！！！！！ // 因此在输出浮点数之前，一定要调用次函数进行修正！ inline double FIX(double x) { return (fabs(x) < EPS) ? 0 : x; } ////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////// //二维矢量运算 bool operator==(Point p1, Point p2) { return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) ); } bool operator!=(Point p1, Point p2) { return ( NEQ(p1.x, p2.x) || NEQ(p1.y, p2.y) ); } bool operator<(Point p1, Point p2) { if (NEQ(p1.x, p2.x)) { return (p1.x < p2.x); } else { return (p1.y < p2.y); } } Point operator+(Point p1, Point p2) { return Point(p1.x + p2.x, p1.y + p2.y); } Point operator-(Point p1, Point p2) { return Point(p1.x - p2.x, p1.y - p2.y); } double operator*(Point p1, Point p2) // 计算叉乘 p1 × p2 { return (p1.x * p2.y - p2.x * p1.y); } double operator&(Point p1, Point p2) { // 计算点积 p1·p2 return (p1.x * p2.x + p1.y * p2.y); } double Norm(Point p) // 计算矢量p的模 { return sqrt(p.x * p.x + p.y * p.y); } // 把矢量p旋转角度angle (弧度表示) // angle > 0表示逆时针旋转 // angle < 0表示顺时针旋转 Point Rotate(Point p, double angle) { Point result; result.x = p.x * cos(angle) - p.y * sin(angle); result.y = p.x * sin(angle) + p.y * cos(angle); return result; } ////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// //三维矢量运算 bool operator==(Point3D p1, Point3D p2) { return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) && EQ(p1.z, p2.z) ); } bool operator<(Point3D p1, Point3D p2) { if (NEQ(p1.x, p2.x)) { return (p1.x < p2.x); } else if (NEQ(p1.y, p2.y)) { return (p1.y < p2.y); } else { return (p1.z < p2.z); } } Point3D operator+(Point3D p1, Point3D p2) { return Point3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z); } Point3D operator-(Point3D p1, Point3D p2) { return Point3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z); } Point3D operator*(Point3D p1, Point3D p2) // 计算叉乘 p1 x p2 { return Point3D(p1.y * p2.z - p1.z * p2.y, p1.z * p2.x - p1.x * p2.z, p1.x * p2.y - p1.y * p2.x ); } double operator&(Point3D p1, Point3D p2) { // 计算点积 p1·p2 return (p1.x * p2.x + p1.y * p2.y + p1.z * p2.z); } double Norm(Point3D p) // 计算矢量p的模 { return sqrt(p.x * p.x + p.y * p.y + p.z * p.z); } ////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////// //点.线段.直线问题 // double Distance(Point p1, Point p2) //2点间的距离 { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } double Distance(Point3D p1, Point3D p2) //2点间的距离,三维 { 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)); } double Distance(Point p, Line L) // 求二维平面上点到直线的距离 { return ( fabs((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) ); } double Distance(Point3D p, Line3D L)// 求三维空间中点到直线的距离 { return ( Norm((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) ); } bool OnLine(Point p, Line L) // 判断二维平面上点p是否在直线L上 { return ZERO( (p - L.p1) * (L.p2 - L.p1) ); } bool OnLine(Point3D p, Line3D L) // 判断三维空间中点p是否在直线L上 { return ZERO( (p - L.p1) * (L.p2 - L.p1) ); } int Relation(Point p, Line L) // 计算点p与直线L的相对关系 ,返回ONLINE,LEFT,RIGHT { double res = (L.p2 - L.p1) * (p - L.p1); if (EQ(res, 0)) { return ONLINE; } else if (res > 0) { return LEFT; } else { return RIGHT; } } bool SameSide(Point p1, Point p2, Line L) // 判断点p1, p2是否在直线L的同侧 { double m1 = (p1 - L.p1) * (L.p2 - L.p1); double m2 = (p2 - L.p1) * (L.p2 - L.p1); return GT(m1 * m2, 0); } bool OnLineSeg(Point p, Line L) // 判断二维平面上点p是否在线段l上 { return ( ZERO( (L.p1 - p) * (L.p2 - p) ) && LEQ((p.x - L.p1.x)*(p.x - L.p2.x), 0) && LEQ((p.y - L.p1.y)*(p.y - L.p2.y), 0) ); } bool OnLineSeg(Point3D p, Line3D L) // 判断三维空间中点p是否在线段l上 { return ( ZERO((L.p1 - p) * (L.p2 - p)) && EQ( Norm(p - L.p1) + Norm(p - L.p2), Norm(L.p2 - L.p1)) ); } Point SymPoint(Point p, Line L) // 求二维平面上点p关于直线L的对称点 { Point result; double a = L.p2.x - L.p1.x; double b = L.p2.y - L.p1.y; double t = ( (p.x - L.p1.x) * a + (p.y - L.p1.y) * b ) / (a*a + b*b); result.x = 2 * L.p1.x + 2 * a * t - p.x; result.y = 2 * L.p1.y + 2 * b * t - p.y; return result; } bool Coplanar(Points3D points) // 判断一个点集中的点是否全部共面 { int i; Point3D p; if (points.size() < 4) return true; p = (points[2] - points[0]) * (points[1] - points[0]); for (i = 3; i < points.size(); i++) { if (! ZERO(p & points[i]) ) return false; } return true; } bool LineIntersect(Line L1, Line L2) // 判断二维的两直线是否相交 { return (! ZERO((L1.p1 - L1.p2)*(L2.p1 - L2.p2)) ); // 是否平行 } bool LineIntersect(Line3D L1, Line3D L2) // 判断三维的两直线是否相交 { Point3D p1 = L1.p1 - L1.p2; Point3D p2 = L2.p1 - L2.p2; Point3D p = p1 * p2; if (ZERO(p)) return false; // 是否平行 p = (L2.p1 - L1.p2) * (L1.p1 - L1.p2); return ZERO(p & L2.p2); // 是否共面 } bool LineSegIntersect(Line L1, Line L2) // 判断二维的两条线段是否相交 { return ( GEQ( max(L1.p1.x, L1.p2.x), min(L2.p1.x, L2.p2.x) ) && GEQ( max(L2.p1.x, L2.p2.x), min(L1.p1.x, L1.p2.x) ) && GEQ( max(L1.p1.y, L1.p2.y), min(L2.p1.y, L2.p2.y) ) && GEQ( max(L2.p1.y, L2.p2.y), min(L1.p1.y, L1.p2.y) ) && LEQ( ((L2.p1 - L1.p1) * (L1.p2 - L1.p1)) * ((L2.p2 - L1.p1) * (L1.p2 - L1.p1)), 0 ) && LEQ( ((L1.p1 - L2.p1) * (L2.p2 - L2.p1)) * ((L1.p2 - L2.p1) * (L2.p2 - L2.p1)), 0 ) ); } bool LineSegIntersect(Line3D L1, Line3D L2) // 判断三维的两条线段是否相交 { // todo return true; } // 计算两条二维直线的交点，结果在参数P中返回 // 返回值说明了两条直线的位置关系: COLINE -- 共线 PARALLEL -- 平行 CROSS -- 相交 int CalCrossPoint(Line L1, Line L2, Point& P) { double A1, B1, C1, A2, B2, C2; A1 = L1.p2.y - L1.p1.y; B1 = L1.p1.x - L1.p2.x; C1 = L1.p2.x * L1.p1.y - L1.p1.x * L1.p2.y; A2 = L2.p2.y - L2.p1.y; B2 = L2.p1.x - L2.p2.x; C2 = L2.p2.x * L2.p1.y - L2.p1.x * L2.p2.y; if (EQ(A1 * B2, B1 * A2)) { if (EQ( (A1 + B1) * C2, (A2 + B2) * C1 )) { return COLINE; } else { return PARALLEL; } } else { P.x = (B2 * C1 - B1 * C2) / (A2 * B1 - A1 * B2); P.y = (A1 * C2 - A2 * C1) / (A2 * B1 - A1 * B2); return CROSS; } } // 计算两条三维直线的交点，结果在参数P中返回 // 返回值说明了两条直线的位置关系 COLINE -- 共线 PARALLEL -- 平行 CROSS -- 相交 NONCOPLANAR -- 不公面 int CalCrossPoint(Line3D L1, Line3D L2, Point3D& P) { // todo return 0; } // 计算点P到直线L的最近点 Point NearestPointToLine(Point P, Line L) { Point result; double a, b, t; a = L.p2.x - L.p1.x; b = L.p2.y - L.p1.y; t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b); result.x = L.p1.x + a * t; result.y = L.p1.y + b * t; return result; } // 计算点P到线段L的最近点 Point NearestPointToLineSeg(Point P, Line L) { Point result; double a, b, t; a = L.p2.x - L.p1.x; b = L.p2.y - L.p1.y; t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b); if ( GEQ(t, 0) && LEQ(t, 1) ) { result.x = L.p1.x + a * t; result.y = L.p1.y + b * t; } else { if ( Norm(P - L.p1) < Norm(P - L.p2) ) { result = L.p1; } else { result = L.p2; } } return result; } // 计算险段L1到线段L2的最短距离 double MinDistance(Line L1, Line L2) { double d1, d2, d3, d4; if (LineSegIntersect(L1, L2)) { return 0; } else { d1 = Norm( NearestPointToLineSeg(L1.p1, L2) - L1.p1 ); d2 = Norm( NearestPointToLineSeg(L1.p2, L2) - L1.p2 ); d3 = Norm( NearestPointToLineSeg(L2.p1, L1) - L2.p1 ); d4 = Norm( NearestPointToLineSeg(L2.p2, L1) - L2.p2 ); return min( min(d1, d2), min(d3, d4) ); } } // 求二维两直线的夹角， // 返回值是0~Pi之间的弧度 double Inclination(Line L1, Line L2) { Point u = L1.p2 - L1.p1; Point v = L2.p2 - L2.p1; return acos( (u & v) / (Norm(u)*Norm(v)) ); } // 求三维两直线的夹角， // 返回值是0~Pi之间的弧度 double Inclination(Line3D L1, Line3D L2) { Point3D u = L1.p2 - L1.p1; Point3D v = L2.p2 - L2.p1; return acos( (u & v) / (Norm(u)*Norm(v)) ); } ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// // 判断两个矩形是否相交 // 如果相邻不算相交 bool Intersect(Rect_2 r1, Rect_2 r2) { return ( max(r1.xl, r2.xl) < min(r1.xh, r2.xh) && max(r1.yl, r2.yl) < min(r1.yh, r2.yh) ); } // 判断矩形r2是否可以放置在矩形r1内 // r2可以任意地旋转 //发现原来的给出的方法过不了OJ上的无归之室这题， //所以用了自己的代码 bool IsContain(Rect r1, Rect r2) //矩形的w>h { if(r1.w >r2.w && r1.h > r2.h) return true; else { double r = sqrt(r2.w*r2.w + r2.h*r2.h) / 2.0; double alpha = atan2(r2.h,r2.w); double sita = asin((r1.h/2.0)/r); double x = r * cos(sita - 2*alpha); double y = r * sin(sita - 2*alpha); if(x < r1.w/2.0 && y < r1.h/2.0 && x > 0 && y > -r1.h/2.0) return true; else return false; } } //////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////// //圆 Point Center(const Circle & C) //圆心 { return C.c; } double Area(const Circle &C) { return pi*C.r*C.r; } double CommonArea(const Circle & A, const Circle & B) //两个圆的公共面积 { double area = 0.0; const Circle & M = (A.r > B.r) ? A : B; const Circle & N = (A.r > B.r) ? B : A; double D = Distance(Center(M), Center(N)); if ((D < M.r + N.r) && (D > M.r - N.r)) { double cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D); double cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D); double alpha = 2.0 * acos(cosM); double beta = 2.0 * acos(cosN); double TM = 0.5 * M.r * M.r * sin(alpha); double TN = 0.5 * N.r * N.r * sin(beta); double FM = (alpha / (2*pi)) * Area(M); double FN = (beta / (2*pi)) * Area(N); area = FM + FN - TM - TN; } else if (D <= M.r - N.r) { area = Area(N); } return area; } bool IsInCircle(const Circle & C, const Rect_2 & rect)//判断圆是否在矩形内(不允许相切) { return (GT(C.c.x - C.r, rect.xl) && LT(C.c.x + C.r, rect.xh) && GT(C.c.y - C.r, rect.yl) && LT(C.c.y + C.r, rect.yh)); } //判断2圆的位置关系 //返回: //BAOHAN = 1; // 大圆包含小圆 //NEIQIE = 2; // 内切 //XIANJIAO = 3; // 相交 //WAIQIE = 4; // 外切 //XIANLI = 5; // 相离 int CirCir(const Circle &c1, const Circle &c2)//判断2圆的位置关系 { double dis = Distance(c1.c,c2.c); if(LT(dis,fabs(c1.r-c2.r))) return BAOHAN; if(EQ(dis,fabs(c1.r-c2.r))) return NEIQIE; if(LT(dis,c1.r+c2.r) && GT(dis,fabs(c1.r-c2.r))) return XIANJIAO; if(EQ(dis,c1.r+c2.r)) return WAIQIE; return XIANLI; } ////////////////////////////////////////////////////////////////////////
1.14结构体表示几何图形
//计算几何(二维) #include <cmath> #include <cstdio> #include <algorithm> using namespace std; typedef double TYPE; #define Abs(x) (((x)>0)?(x):(-(x))) #define Sgn(x) (((x)<0)?(-1):(1)) #define Max(a,b) (((a)>(b))?(a):(b)) #define Min(a,b) (((a)<(b))?(a):(b)) #define Epsilon 1e-8 #define Infinity 1e+10 #define PI acos(-1.0)//3.14159265358979323846 TYPE Deg2Rad(TYPE deg){return (deg * PI / 180.0);} TYPE Rad2Deg(TYPE rad){return (rad * 180.0 / PI);} TYPE Sin(TYPE deg){return sin(Deg2Rad(deg));} TYPE Cos(TYPE deg){return cos(Deg2Rad(deg));} TYPE ArcSin(TYPE val){return Rad2Deg(asin(val));} TYPE ArcCos(TYPE val){return Rad2Deg(acos(val));} TYPE Sqrt(TYPE val){return sqrt(val);} //点 struct POINT { TYPE x; TYPE y; POINT() : x(0), y(0) {}; POINT(TYPE _x_, TYPE _y_) : x(_x_), y(_y_) {}; }; // 两个点的距离 TYPE Distance(const POINT & a, const POINT & b) { return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); } //线段 struct SEG { POINT a; //起点 POINT b; //终点 SEG() {}; SEG(POINT _a_, POINT _b_):a(_a_),b(_b_) {}; }; //直线(两点式) struct LINE { POINT a; POINT b; LINE() {}; LINE(POINT _a_, POINT _b_) : a(_a_), b(_b_) {}; }; //直线(一般式) struct LINE2 { TYPE A,B,C; LINE2() {}; LINE2(TYPE _A_, TYPE _B_, TYPE _C_) : A(_A_), B(_B_), C(_C_) {}; }; //两点式化一般式 LINE2 Line2line(const LINE & L) // y=kx+c k=y/x { LINE2 L2; L2.A = L.b.y - L.a.y; L2.B = L.a.x - L.b.x; L2.C = L.b.x * L.a.y - L.a.x * L.b.y; return L2; } // 引用返回直线 Ax + By + C =0 的系数 void Coefficient(const LINE & L, TYPE & A, TYPE & B, TYPE & C) { A = L.b.y - L.a.y; B = L.a.x - L.b.x; C = L.b.x * L.a.y - L.a.x * L.b.y; } void Coefficient(const POINT & p,const TYPE a,TYPE & A,TYPE & B,TYPE & C) { A = Cos(a); B = Sin(a); C = - (p.y * B + p.x * A); } //判等(值，点，直线) bool IsEqual(TYPE a, TYPE b) { return (Abs(a - b) <Epsilon); } bool IsEqual(const POINT & a, const POINT & b) { return (IsEqual(a.x, b.x) && IsEqual(a.y, b.y)); } bool IsEqual(const LINE & A, const LINE & B) { TYPE A1, B1, C1; TYPE A2, B2, C2; Coefficient(A, A1, B1, C1); Coefficient(B, A2, B2, C2); return IsEqual(A1 * B2, A2 * B1) && IsEqual(A1 * C2, A2 * C1) && IsEqual(B1 * C2, B2 * C1); } // 矩形 struct RECT { POINT a; // 左下点 POINT b; // 右上点 RECT() {}; RECT(const POINT & _a_, const POINT & _b_) { a = _a_; b = _b_; } }; //矩形化标准 RECT Stdrect(const RECT & q) { TYPE t; RECT p=q; if(p.a.x > p.b.x) swap(p.a.x , p.b.x); if(p.a.y > p.b.y) swap(p.a.y , p.b.y); return p; } //根据下标返回矩形的边 SEG Edge(const RECT & rect, int idx) { SEG edge; while (idx < 0) idx += 4; switch (idx % 4) { case 0: //下边 edge.a = rect.a; edge.b = POINT(rect.b.x, rect.a.y); break; case 1: //右边 edge.a = POINT(rect.b.x, rect.a.y); edge.b = rect.b; break; case 2: //上边 edge.a = rect.b; edge.b = POINT(rect.a.x, rect.b.y); break; case 3: //左边 edge.a = POINT(rect.a.x, rect.b.y); edge.b = rect.a; break; default: break; } return edge; } //矩形的面积 TYPE Area(const RECT & rect) { return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y); } //两个矩形的公共面积 TYPE CommonArea(const RECT & A, const RECT & B) { TYPE area = 0.0; POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y)); POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y)); if( (LL.x <= UR.x) && (LL.y <= UR.y) ) { area = Area(RECT(LL, UR)); } return area; } //判断圆是否在矩形内(不允许相切) bool IsInCircle(const CIRCLE & circle, const RECT & rect) { return (circle.x - circle.r > rect.a.x) && (circle.x + circle.r < rect.b.x) && (circle.y - circle.r > rect.a.y) && (circle.y + circle.r < rect.b.y); } //判断矩形是否在圆内(不允许相切) bool IsInRect(const CIRCLE & circle, const RECT & rect) { POINT c,d; c.x=rect.a.x; c.y=rect.b.y; d.x=rect.b.x; d.y=rect.a.y; return (Distance( Center(circle) , rect.a ) < circle.r) && (Distance( Center(circle) , rect.b ) < circle.r) && (Distance( Center(circle) , c ) < circle.r) && (Distance( Center(circle) , d ) < circle.r); } //判断矩形是否与圆相离(不允许相切) bool Isoutside(const CIRCLE & circle, const RECT & rect) { POINT c,d; c.x=rect.a.x; c.y=rect.b.y; d.x=rect.b.x; d.y=rect.a.y; return (Distance( Center(circle) , rect.a ) > circle.r) && (Distance( Center(circle) , rect.b ) > circle.r) && (Distance( Center(circle) , c ) > circle.r) && (Distance( Center(circle) , d ) > circle.r) && (rect.a.x > circle.x || circle.x > rect.b.x || rect.a.y > circle.y || circle.y > rect.b.y) || ((circle.x - circle.r > rect.b.x) || (circle.x + circle.r < rect.a.x) || (circle.y - circle.r > rect.b.y) || (circle.y + circle.r < rect.a.y)); }
求多边形最大宽度
#include <bits/stdc++.h> using namespace std; #define ll long long #define lowbit(a) ((a) & -(a)) #define clean(a, b) memset(a, b, sizeof(a)) const int mod = 1e9 + 7; const double inf = 1e50; const int maxn = 1e5 + 9; ////////////////////////////////////////////////////////////////////////// int n; struct node { int x, y; node() {} node(int _x, int _y){ x = _x, y = _y;} node operator-(node b) { return node(x - b.x, y - b.y); } double len(){ return hypot(x, y); } } a[maxn]; double cross(node a, node b) { return a.x * b.y - a.y * b.x; } double dist(node p, node a, node b) { return cross(p - a, b - a) / (b - a).len(); } int main() { scanf("%d", &n); for (int i = 1; i <= n; i++) scanf("%d %d", &a[i].x, &a[i].y); double ans = inf; for (int i = 1; i <= n; i++) { for (int j = 1; j < i; j++) { double l = inf, r = -inf; for (int k = 1; k <= n; k++) { l = min(l, dist(a[k], a[i], a[j])); r = max(r, dist(a[k], a[i], a[j])); } r -= l; ans = min(ans, r); } } printf("%.10f\n", ans); return 0; }
浙公网安备 33010602011771号