陈家泽 计算几何模板
7 计 算 几 何
7.1 几 何 基 础
7.1.1 求 两 向 量 的 叉积
int cross(int x1,int y1,int x2,int y2) //int型 { return x1*y2-x2*y1; } double cross(double x1,double y1,double x2,double y2) //double型 { return x1*y2-x2*y1; } double cross(Point a,Point b) // 向量叉积 { return a.x*b.y-b.x*a.y; }
7.1.2 求 平 面 两 点 欧 氏 距离
double distance(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); }
7.1.5 判 断double 型 变 量 的 符号（ps：重点，替换kuangbin sgn函数）
const double eps=1e-8; int dlcmp(double x) { return x<-eps?-1:x>eps; }
7.1.7 求 两 线 段 间 最 短 距离
double dis_segments(Segment seg1,Segment seg2) { double m1=dis_point_segment(seg1.s,seg2); double m2=dis_point_segment(seg1.e,seg2); double m3=dis_point_segment(seg2.s,seg1); double m4=dis_point_segment(seg2.e,seg1); return min(min(m1,m2),min(m3,m4)); }
7.1.8 判 断 两 直 线 是 否 相 交 （ 共 线 及 平 行 ， 若 相 交 则 求 出 交 点 ） （请认真考虑一下直线与线段，以及与kuangbin模板里的&操作符）
#include <math.h> #define eps 1e-8 struct Point { double x; double y; }; struct Line { Point s; Point e; Point v; }; Line l1,l2; double x,y交点坐标;// int line_intersect(Line l1,Line l2) { Point vec; double r; vec.x=l1.s.x-l2.s.x; vec.y=l1.s.y-l2.s.y; if (fabs(l1.v.x*l2.v.y-l2.v.x*l1.v.y)<eps) { if (fabs(l1.v.x*vec.y-vec.x*l1.v.y)<eps) return 2; //共线 else return 0; //平行 } else { r=((l1.s.x-l2.s.x)*l2.v.y-(l1.s.y-l2.s.y)*l2.v.x) /(l1.v.x*l2.v.y-l1.v.y*l2.v.x); x=-r*l1.v.x+l1.s.x; y=-r*l1.v.y+l1.s.y; return 1; //相交 } }
7.2 多 边 形
7.2.1 判 断 线 段 是 否 与 矩 形 相 交 （ 包 括 线 段 在 矩 形 内 部 ）
struct Rectangle//这好像是矩形 { Point lt;//lefttop Point rb;//rightbottom }; int segement_rectangle_intersect(Segment l,Rectangle r) { Segment d1,d2;//retangle ’s diagonal d1.s=r.lt; d1.e=r.rb; d2.s.x=d1.e.x; d2.s.y=d1.s.y; d2.e.x=d1.s.x; d2.e.y=d1.e.y; if (l.s.x>=r.lt.x&&l.s.x<=r.rb.x&& l.s.y<=r.lt.y&&l.s.y>=r.rb.y|| l.e.x>=r.lt.x&&l.e.x<=r.rb.x&& l.e.y<=r.lt.y&&l.e.y>=r.rb.y) return 1; if (segment_intersect(l,d1)|| segment_intersect(l,d2)) //segment_intersect(endpoint inclusive) return 1; return 0; }
7.2.3 求 简 单 多 边 形 重 心
Point get_center(Point pt[],int n) { double sum,area; Point res(0,0),o(0,0); int i; sum=0; for (i=0; i<n; i++) { area=cross(pt[i]-o,pt[(i+1)%n]-o); res=res+(pt[i]+pt[(i+1)%n])/3*area; sum+=area; } res=res/sum; return res; }
7.2.7 判 断 两 凸 多 边 形 是 否 相 交 （graham scan 求 凸 包+ 枚 举 边 、 点)
#include <iostream > #include <cstdio > #include <algorithm > #define N 110 using namespace std; struct Point { int x,y; }; struct Polygon { Point p[N]; int n; }; Point pt[N]; int stack[N]; int cross(Point a,Point b,Point s) { return (a.x-s.x)*(b.y-s.y)-(b.x-s.x)*(a.y-s.y); } int dist(Point a,Point b) { return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y); } int cmp(Point a,Point b) { if (cross(a,b,pt[0]) >0||cross(a,b,pt[0])==0&&dist(a,pt[0])<dist(b,pt[0])) return 1; else return 0; } int on_segment(Point s,Point e,Point o) { if(cross(s,e,o)==0&& o.x>=min(s.x,e.x)&& o.x<=max(s.x,e.x)&& o.y>=min(s.y,e.y)&& o.y<=max(s.y,e.y)) return 1; else return 0; } int graham_scan(int n) { int i,top,t; if (n<=1) { stack[0]=0; return n; } for (t=0,i=1; i<n; i++) if (pt[i].y<pt[t].y||pt[i].y==pt[t].y&&pt[i].x<pt[t].x) t=i; swap(pt[0],pt[t]); sort(pt+1,pt+n,cmp); top=2; for (i=0; i<2; i++) stack[i]=i; for (i=2; i<n; i++) { while (top >1&&cross(pt[stack[top -1]],pt[i],pt[stack[top -2]]) <=0) top--; stack[top++]=i; } return top; } int segment_intersect(Point s1,Point e1,Point s2,Point e2) { if (max(s1.x,e1.x)>=min(s2.x,e2.x)&& max(s2.x,e2.x)>=min(s1.x,e1.x)&& max(s1.y,e1.y)>=min(s2.y,e2.y)&& max(s2.y,e2.y)>=min(s1.y,e2.y)&& (double)cross(s2,e1,s1)*(double)cross(e2,e1,s1)<=0&& (double)cross(s1,e2,s2)*(double)cross(e1,e2,s2)<=0) return 1; else return 0; } int point_inside(Point o,Polygon pln) { int i,a1=0,a2=0,n=pln.n; pln.p[n]=pln.p[0]; for (i=0; i<n; i++) a1+=abs(cross(pln.p[i],pln.p[i+1],o)); for (i=1; i<n; i++) a2+=abs(cross(pln.p[i],pln.p[i+1],pln.p[0])); if (a1==a2) return 1; else return 0; } int convex_polygon_intersect(Polygon pln1,Polygon pln2) { int i,j; pln1.p[pln1.n]=pln1.p[0]; pln2.p[pln2.n]=pln2.p[0]; for (i=0; i<pln1.n; i++) for (j=0; j<pln2.n; j++) if (segment_intersect(pln1.p[i],pln1.p[i+1],pln2.p[j],pln2.p[j+1])) return 1; if (point_inside(pln1.p[0],pln2)||point_inside(pln2.p[0],pln1)) return 1; return 0; } int main() { int n,m,i,vertexnum,ans; Polygon pln1,pln2; while (cin>>n>>m,n||m) { for (i=0; i<n; i++) cin>>pt[i].x>>pt[i].y; vertexnum=graham_scan(n); pln1.n=vertexnum; for (i=0; i<vertexnum; i++) pln1.p[i]=pt[stack[i]]; for (i=0; i<m; i++) cin>>pt[i].x>>pt[i].y; vertexnum=graham_scan(m); pln2.n=vertexnum; for (i=0; i<vertexnum; i++) pln2.p[i]=pt[stack[i]]; if (pln1.n==1&&pln2.n==1) ans=1; else if (pln1.n==1&&pln2.n==2) { if (on_segment(pln2.p[0],pln2.p[1],pln1.p[0])) ans=0; else ans=1; } else if (pln1.n==2&&pln2.n==1) { if (on_segment(pln1.p[0],pln1.p[1],pln2.p[0])) ans=0; else ans=1; } else if (pln1.n==2&&pln2.n==2) { if (segment_intersect(pln1.p[0],pln1.p[1],pln2.p[0],pln2.p[1])) ans=0; else ans=1; } else if (pln1.n==1) { if (point_inside(pln1.p[0],pln2)) ans=0; else ans=1; } else if (pln2.n==1) { if (point_inside(pln2.p[0],pln1)) ans=0; else ans=1; } else { if (convex_polygon_intersect(pln1,pln2)==0) ans=1; else ans=0; } if (ans==0) cout <<"NO"<<endl; else cout <<"YES"<<endl; } return 0; }
7.3.3 两 个 简 单 多 边 形 求 面 积 并 、交
#include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #define eps 1e-8 #define N 550 using namespace std; struct Point { double x,y; Point () {} Point (double a,double b):x(a),y(b) {} Point operator - (const Point a) const { return Point(x-a.x,y-a.y); } }; Point zero(0,0); int dlcmp(double x) { return x<-eps?-1:x>eps; } double cross(Point v1,Point v2) { return v1.x*v2.y-v2.x*v1.y; } struct Polygon { Point p[N]; int n; Polygon():n(0) {} void clear() { n=0; } void add(Point a) { p[n++]=a; } double area() { double res=0; for (int i=1; i<n-1; i++) res+=cross(p[i]-p[0],p[i+1]-p[0]); return fabs(res/2); } }; Polygon A,B,rec; Point line_intersect(Point s1,Point e1,Point s2,Point e2) //两直线交点 { Point v1,v2,res; double r; v1=s1-e1; v2=s2-e2; r=((s1.x-s2.x)*v2.y-(s1.y-s2.y)*v2.x)/(v1.x*v2.y-v1.y*v2.x); res.x=-r*v1.x+s1.x; res.y=-r*v1.y+s1.y; return res; } void cut(Point s,Point e) //半平面交 { int i,j,d1,d2; Polygon ker; for (i=0; i<rec.n; i++) { j=(i+1)%rec.n; d1=dlcmp(cross(e-s,rec.p[i]-s)); d2=dlcmp(cross(e-s,rec.p[j]-s)); if (d1>=0) ker.add(rec.p[i]); if (d1*d2<0) ker.add(line_intersect(s,e,rec.p[i],rec.p[j])); } rec=ker; } double calc(Point p1,Point p2,Point q1,Point q2) { int dp=dlcmp(cross(p1,p2)),dq=dlcmp(cross(q1,q2)); int sgn=dp*dq; if (sgn==0) return 0; rec.clear(); rec.add(zero); rec.add(p1); rec.add(p2); if (dp<0) swap(rec.p[1],rec.p[2]); if (dq>0) { cut(zero,q1); cut(q1,q2); cut(q2,zero); } else { cut(zero,q2); cut(q2,q1); cut(q1,zero); } return sgn*rec.area(); } double solve() { double res=A.area()+B.area(); double sum=0; int i,j; //对两个多边形三角剖分，分别求两个三角形的面积交 for (i=0; i<A.n; i++) for (j=0; j<B.n; j++) sum+=calc(A.p[i],A.p[(i+1)%A.n],B.p[j],B.p[(j+1)%B.n]); res-=fabs(sum); //fabs(sum)为两个多边形的面积交 return res; //面积并 } //hdu3060(题目数据有误) int main() { int n,m,i; Point pt; double ans; while (scanf("%d%d",&n,&m)!=EOF) { A.clear(); B.clear(); for (i=0; i<n; i++) { scanf("%lf%lf",&pt.x,&pt.y); A.add(pt); } for (i=0; i<m; i++) { scanf("%lf%lf",&pt.x,&pt.y); B.add(pt); } ans=solve(); printf("%.2f\n",ans+eps); } return 0; }
7.4 圆
7.4.1 点类
struct Point { double x,y; Point() {} Point(double a,double b):x(a),y(b) {} Point operator + (const Point a) const { return Point(x+a.x,y+a.y); } Point operator - (const Point a) const { return Point(x-a.x,y-a.y); } Point operator * (const double a) const { return Point(x*a,y*a); } Point operator / (const double a) const { return Point(x/a,y/a); } bool operator < (const Point a) const { if (dlcmp(x-a.x)==0) return dlcmp(x-a.y)<0; else return dlcmp(x-a.x)<0; } bool operator == (const Point a) const { return !dlcmp(x-a.x)&&!dlcmp(y-a.y); } //向量长度定为d Point trunc(double d) { double dis(Point,Point); double len=dis(*this,Point(0,0)); return Point(x*d/len,y*d/len); } //坐标逆时针旋转a度 Point rotate(double a) { return Point(x*cos(a)-y*sin(a),y*cos(a)+x*sin(a)); } }; double dis(Point a,Point b) { return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)); } double cross(Point a,Point b,Point s) { double x1=a.x-s.x,y1=a.y-s.y; double x2=b.x-s.x,y2=b.y-s.y; return x1*y2-x2*y1; } double cross(Point a,Point b) { return a.x*b.y-b.x*a.y; } double dot(Point a,Point b,Point s) { double x1=a.x-s.x,y1=a.y-s.y; double x2=b.x-s.x,y2=b.y-s.y; return x1*x2+y1*y2; } double dot(Point a,Point b) { return a.x*b.x+a.y*b.y; }
7.4.2 圆 类
struct Circle { Point o; double r; Circle() {} Circle(Point a,double l):o(a),r(l) {} double area() { return sqr(r)*PI; } }; //判断圆a是否含于圆b int inner_circle(Circle a,Circle b) { if (dlcmp(a.r-b.r)>0) return 0; return dlcmp(dis(a.o,b.o)+a.r-b.r)<=0; } //以base点为基点，极角排序，排序前base需赋初值 Point base; int cmp(const Point a,const Point b) { return atan2(a.y-base.y,a.x-base.x)<atan2(b.y-base.y,b.x-base.x); } //向量a,b的夹角 double vec_angle(Point a,Point b) { double tmp=dot(a,b)/(dis(a,Point(0,0))*dis(b,Point(0,0))); if (dlcmp(tmp -1) >=0) tmp=1; if (dlcmp(tmp+1) <=0) tmp=-1; return acos(tmp); } //计算由a到b逆时针方向的弓形面积 double arc_area(Point a,Point b,Circle c) { double theta=vec_angle(a-c.o,b-c.o); double sf=sqr(c.r)*theta/2.0; double st=sqr(c.r)*sin(theta)/2.0; if (dlcmp(cross(a,b,c.o))>0) return sf-st; else return c.area()-sf+st; } double arc_area(double th,double r) { return 0.5*sqr(r)*(th-sin(th)); }
7.4.3 圆 面 积 交 、并
//求两圆交点，排除相切的情况，不考虑内含 int inter_circle_or(Circle c1,Circle c2,Point &p1,Point &p2) { double len=dis(c1.o,c2.o); if (dlcmp(len-c1.r-c2.r)>=0) return 0; double s=(sqr(c1.r)-sqr(c2.r)+sqr(len))/len/2; double h=sqrt(sqr(c1.r)-sqr(s)); Point vec=c2.o-c1.o; Point p0=c1.o+vec.trunc(s); //圆类里的函数，可能和kuangbin相同 p1=p0+vec.rotate(PI/2).trunc(h); p2=p0-vec.rotate(PI/2).trunc(h); return 1; }
7.4.4 简 单 多 边 形 与 圆 求 面 积 交
#include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #define N 200 #define eps 1e-8 using namespace std; const double PI=acos(-1.0); struct Point { double x,y; }; Point pt[N]; int n; int dlcmp(double x) { return x<-eps?-1:x>eps; } double sqr(double x) { return x*x; } double dis(Point a,Point b) { return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)); } double outer(Point a,Point b,Point c) { return (a.x-c.x)*(b.y-c.y)-(a.y-c.y)*(b.x-c.x); } double inner(Point a,Point b,Point c) { return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y); } double calc_area(Point a,Point b,Point c,double r) { double A,B,C,x,y,tS; A=dis(b,c); B=dis(a,c); C=dis(b,a); if (A<r&&B<r) return outer(a,b,c)/2; else if (A<r&&B>=r) { x=(inner(a,c,b)+sqrt(sqr(r)*sqr(C)-sqr(outer(a,c,b))))/C; tS=outer(a,b,c)/2; return asin(tS*(1-x/C)*2/r/B)*sqr(r)/2+tS*x/C; } else if (A>=r&&B<r) { y=(inner(b,c,a)+sqrt(sqr(r)*sqr(C)-sqr(outer(b,c,a))))/C; tS=outer(a,b,c)/2; return asin(tS*(1-y/C)*2/r/A)*sqr(r)/2+tS*y/C; } else if (fabs(outer(a,b,c))>=r*C||inner(b,c,a) <=0||inner(a,c,b)<=0) { if (inner(a,b,c)<0) { if (outer(a,b,c)<0) return (-PI-asin(outer(a,b,c)/A/B))*sqr(r)/2; else return (PI-asin(outer(a,b,c)/A/B))*sqr(r)/2; } else return asin(outer(a,b,c)/A/B)*sqr(r)/2; } else { x=(inner(a,c,b)+sqrt(sqr(r)*sqr(C)-sqr(outer(a,c,b))))/C; y=(inner(b,c,a)+sqrt(sqr(r)*sqr(C)-sqr(outer(b,c,a))))/C; tS=outer(a,b,c)/2; return (asin(tS*(1-x/C)*2/r/B)+asin(tS*(1-y/C)*2/r/A))*sqr(r)/2+tS*((y+x)/C-1); } } //计算一般多边形与圆的交面积（将多边形划分为三角形，然后有向三角形与圆求有向面积交） double solve(Point o,double r) { int i,j; double res,sum; Point tri[3]; res=0; for (i=1;i<n-1;i++) { tri[0]=pt[0]; tri[1]=pt[i]; tri[2]=pt[i+1]; sum=0; for (j=0;j<3;j++) sum+=calc_area(tri[j],tri[(j+1)%3],o,r); //sum为三角形与圆交的有向面积 res+=sum; } return fabs(res); } //poj3675 int main() { double x0,y0,v,vx,vy,g,r,t,theta ,ans; Point o; int i; o.x=o.y=0; while (scanf("%lf",&r)!=EOF) { scanf("%d",&n); for (i=0;i<n;i++) scanf("%lf%lf",&pt[i].x,&pt[i].y); ans=solve(o,r); printf("%.2f\n",ans); } return 0; }
7.4.5 求 线 段 与 圆 的 交 点
若求直线与圆的交点类似，无需讨论k1、k2的取值范围
int inter_circle_segment(Circle c,Point a,Point b,Point &p1,Point &p2) { Point vec=b-a; double A=sqr(vec.x)+sqr(vec.y); double B=2*(vec.x*(a.x-c.o.x)+vec.y*(a.y-c.o.y)); double C=sqr(a.x-c.o.x)+sqr(a.y-c.o.y)-sqr(c.r); double delta=sqr(B)-4*A*C; if (dlcmp(delta)<0) return 0; double k1=(-B-sqrt(fabs(delta)))/(2*A); double k2=(-B+sqrt(fabs(delta)))/(2*A); int res=0; if (dlcmp(k1) >=0&&dlcmp(k1-1) <=0) { res++; p1=a+vec*k1; } if (dlcmp(k2) >=0&&dlcmp(k2-1) <=0) { res++; if (res==1) p1=a+vec*k2; else p2=a+vec*k2; } return res; }
7.4.6 求两圆公切线
//求两相离的圆的两条内共切线 void get_InCommonTangent(Circle c1,Circle c2,Point &s1,Point &e1,Point &s2,Point &e2) { double l=dis(c1.o,c2.o); double d=l*c1.r/(c1.r+c2.r); double tmp=c1.r/d; tmp=fix(tmp); double theta=acos(tmp); Point vec=c2.o-c1.o; vec=vec.trunc(c1.r); s1=c1.o+vec.rotate(theta); s2=c1.o+vec.rotate(-theta); vec=c1.o-c2.o; vec=vec.trunc(c2.r); e1=c2.o+vec.rotate(theta); e2=c2.o+vec.rotate(-theta); } //求两相离的圆的两条外公切线 void get_OutCommonTangent(Circle c1,Circle c2,Point &s1,Point &e1,Point &s2,Point &e2) { double l=dis(c1.o,c2.o); double d=fabs(c1.r-c2.r); double theta=acos(d/l); if (dlcmp(c1.r-c2.r)>0) swap(c1,c2); Point vec=c1.o-c2.o; vec=vec.trunc(c1.r); s1=c1.o+vec.rotate(theta); s2=c1.o+vec.rotate(-theta); vec=vec.trunc(c2.r); e1=c2.o+vec.rotate(theta); e2=c2.o+vec.rotate(-theta); }
7.4.7 最 小 圆 覆盖（平面上n个点，求一个半径最小的圆，能够覆盖所有的点。）
#include <Bits/stdc++.h> using namespace std; #define eps 1e-8 #define MAX_P 2000 struct Point { double x,y; Point operator -(Point &a) { Point t; t.x=x-a.x; t.y=y-a.y; return t; } }; struct Circle { double r; Point center; }; struct Triangle { Point t[3]; }; Point pt[MAX_P]; //点集 Circle c; //最小圆 double distance(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } double cross(Point a,Point b) { return a.x*b.y-b.x*a.y; } double triangle_area(Triangle tri) //三角形距离 { Point v1=tri.t[1]-tri.t[0]; Point v2=tri.t[2]-tri.t[0]; return fabs(cross(v1,v2))/2; } Circle circumcircle_triangle(Triangle tri) //三角形外接圆 { Circle res; double a,b,c,c1,c2; double xA,yA,xB,yB,xC,yC; a=distance(tri.t[0],tri.t[1]); b=distance(tri.t[1],tri.t[2]); c=distance(tri.t[2],tri.t[0]); res.r=a*b*c/triangle_area(tri)/4; xA=tri.t[0].x; yA=tri.t[0].y; xB=tri.t[1].x; yB=tri.t[1].y; xC=tri.t[2].x; yC=tri.t[2].y; c1=(xA*xA+yA*yA-xB*xB-yB*yB)/2; c2 = (xA*xA+yA*yA-xC*xC-yC*yC)/2; res.center.x=(c1*(yA-yC)-c2*(yA-yB))/ ((xA-xB)*(yA-yC)-(xA-xC)*(yA-yB)); res.center.y = (c1*(xA-xC)-c2*(xA-xB))/ ((yA-yB)*(xA-xC)-(yA-yC)*(xA-xB)); return res; } Circle mincircle_triangle(int trinum,Triangle tri) { Circle res; if (trinum==0) res.r=-2; else if (trinum==1) { res.center=tri.t[0]; res.r=0; } else if (trinum==2) { res.center.x=(tri.t[0].x+tri.t[1].x)/2; res.center.y=(tri.t[0].y+tri.t[1].y)/2; res.r=distance(tri.t[0],tri.t[1])/2; } else if (trinum==3) res=circumcircle_triangle(tri); return res; } void mincircle_pointset(int m,int trinum,Triangle tri) //求点集的最小覆盖圆 { int i,j; Point tmp; c=mincircle_triangle(trinum,tri); if (trinum==3) return; for (i=0; i<m; i++) if (distance(pt[i],c.center)>c.r) { tri.t[trinum]=pt[i]; mincircle_pointset(i,trinum+1,tri); tmp=pt[i]; for (j=i; j>=1; j--) pt[j]=pt[j-1]; pt[0]=tmp; } } int main() { int n,i,f1,f2; Triangle tri; while (scanf("%d%d%d",&f1,&f2,&n)!=EOF) { for (i=0; i<n; i++) scanf("%lf%lf",&pt[i].x,&pt[i].y); mincircle_pointset(n,0,tri); printf("%lf %lf %lf\n",c.center.x,c.center.y,c.r); } return 0; }
7.4.8 单 位 圆 覆盖（一个单位圆最多能覆盖平面上多少点）
#include <math.h> #include <bits/stdc++.h> using namespace std; #define eps 1e-8 #define MAX_P 505 const double r=1.0;//单位圆半径 struct Point { double x,y; }; Point pt[MAX_P]; double distance(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); //sqrt函数速度较慢，应尽量避免出现，此处可优化为距离的平方和的形式 } Point find_center(Point a,Point b) { Point v,mid,center; double d,s,ang; v.x=a.x-b.x; v.y=a.y-b.y; mid.x=(a.x+b.x)/2; mid.y=(a.y+b.y)/2; d=distance(a,mid); s=sqrt(r*r-d*d); //优化为s=sqrt(r*r-d); if (fabs(v.y)<eps) { center.x=mid.x; center.y=mid.y+s; } else { ang=atan(-v.x/v.y); center.x=mid.x+s*cos(ang); center.y=mid.y+s*sin(ang); } return center; } int main() { int n,i,j,k,ans,cnt; double tmp; Point center; while (~scanf("%d",&n),n) { for (i=0; i<n; i++) scanf("%lf%lf",&pt[i].x,&pt[i].y); ans=1; for (i=0; i<n; i++) for (j=i+1; j<n; j++) { if (distance(pt[i],pt[j])>2*r) //优化为distance(pt[i],pt[j])>2*2*r*r continue; cnt=0; center=find_center(pt[i],pt[j]); for (k=0; k<n; k++) if (distance(pt[k],center)<=r+eps) cnt++; if (ans<cnt) ans=cnt; } printf("%d\n",ans); } return 0; }
7.5 模 拟 退 火
7.5.1 求 多 边 形 费 马点
#include <iostream > #include <cstdio > #include <cmath > #define eps 1e-6 #define N 105 using namespace std; struct Point { double x,y; }; double point_dis(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } double sum_dis(Point pt[],int n,Point o) { double res=0; int i; for (i=0;i<n;i++) res+=point_dis(pt[i],o); return res; } double polygon_Fermatpoint(Point pln[],int n) { Point cp,np,tmp; double min,step,d; int flag; cp=pln[0]; //cp保存当前更新后最优的费马点 min=sum_dis(pln,n,cp); step=10000; //选取坐标范围的最大值 while (step>eps) { flag=1; while (flag) { flag=0; np=cp; tmp=cp,tmp.x+=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp, flag=1; tmp=cp,tmp.x-=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; tmp=cp,tmp.y+=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; tmp=cp,tmp.y-=step; d=sum_dis(pln,n,tmp); if (min>d) min=d, np=tmp,flag=1; cp=np; } step*=0.98; //系数根据精度要求修改 } return min; } int main() { int n,i; double min; Point pt[N]; cin>>n; for (i=0;i<n;i++) cin>>pt[i].x>>pt[i].y; min=polygon_Fermatpoint(pt,n); printf("%.0f\n",min); return 0; }
7.5.2 最 小 球 覆盖（求能覆盖所有点的最小球的半径。）
#include <iostream > #include <cstdio > #include <cmath > #define oo 1e20 #define eps 1e-10 #define N 105 using namespace std; struct Point { double x,y,z; }; double dis(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z)); } int max_dis(Point pt[],int n,Point o) { int i,res; double max,tmp; max=0; res=0; for (i=0; i<n; i++) { tmp=dis(pt[i],o); if (max<tmp) { max=tmp; res=i; } } return res; } int main() { Point pt[N],o; int n,i,t; double dx,dy,dz,step,r,tmp; cin>>n; for (i=0; i<n; i++) cin>>pt[i].x>>pt[i].y>>pt[i].z; step=10000; //step选取最大的坐标范围 r=oo; if (n==1) { o.x=pt[0].x; o.y=pt[0].y; o.z=pt[0].z; } else { o.x=o.y=o.z=0; while (step>eps) { t=max_dis(pt,n,o); tmp=dis(pt[t],o); if (r>tmp) r=tmp; dx=(pt[t].x-o.x)/tmp; dy=(pt[t].y-o.y)/tmp; dz=(pt[t].z-o.z)/tmp; o.x+=step*dx; o.y+=step*dy; o.z+=step*dz; step*=0.9993; //系数的选取根据具体精度调整 } } printf("%.6f %.6f %.6f\n",o.x,o.y,o.z); return 0; }
