1.最小的包围圆，将所有的点包围起来。（hdu 3932）最小覆盖圆算法地址：http://soft.cs.tsinghua.edu.cn/blog/?q=node/1066
问题的背景提出：考察固定在工作平台上的一直机械手，要捡起散落在不同位置的多个零件，并送到别的地方。那么，这只机械手的底座应该选在哪里呢？根据直觉，应该选在机械手需够着的那些位置的“中心”。准确地讲，也就是包围这些点的那个最小圆的圆心----该位置的好处是，可使机械手的底座到它需要够着的那些点的最大距离最小化。于是可得如下问题：给定由平面上n个点所组成的一个集合P（对应于机械手需要够着的工作平台的那些位置），试找出P的最小包围圆（smallest enclosing disc）----亦即，包含P中所有点、半径最小的那个圆。这个最小包围圆必然是唯一的。
算法介绍：我们本次算法的设计是基于这样一个简单直观的性质：在既定的给定点条件下，如果引入一张新的半平面，只要此前的最优解顶点（即唯一确定最小包围圆的几个关键顶点）能够包含于其中，则不必对此最优解进行修改，亦即此亦为新点集的最优解；否则，新的最优解顶点必然位于这个新的半空间的边界上。
定理可以通过反证法证明。
于是，基于此性质，我们便可得到一个类似于线性规划算法的随机增量式算法。定义Di为相对于pi的最小包围圆。此算法实现的关键在于对于pi∉Di-1时的处理。显然，如果pi∈Di-1，则Di= Di-1；否则，需要对Di另外更新。而且，Di的组成必然包含了pi；因此，此种情况下的最小包围圆是过pi点且覆盖点集{ p1 ，p2 ，p3 ……pi-1}的最小包围圆。则仿照上述处理的思路，Di={ p1 ，pi }，逐个判断点集{ p2 ，p3 ……pi-1 }，如果存在pj∉ Di，则Di={pj，pi }。同时，再依次对点集{ p1 ，p2 ，p3 ……pj-1 }判断是否满足pk∈Di，若有不满足，则Di={pk ，pj，pi }。由于，三点唯一地确定一个圆，故而，只需在此基础上判断其他的点是否位于此包围圆内，不停地更新pk。当最内层循环完成时，退出循环，转而更新pj；当次内层循环结束时，退出循环，更新pi。当i=n时，表明对所有的顶点均已处理过 ，此时的Dn即表示覆盖了给定n个点的最小包围圆。
问题：找出一个点使得这个店到n个点的最长距离最短，即求最小覆盖圆的半径。
1 #include<cstdio> 2 #include<iostream> 3 #include<cstring> 4 #include<algorithm> 5 #include<cmath> 6 using namespace std; 7 8 const double eps = 1e-8; 9 10 struct node{ 11 double x,y; 12 }; 13 14 struct node p[1005],central; 15 double R; 16 int n; 17 18 double dist(struct node a,struct node b){ 19 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 20 } 21 22 23 //求外接圆圆心,根据三边相等 24 struct node circumcenter(struct node a,struct node b,struct node c) 25 { 26 double a1=b.x-a.x, b1=b.y-a.y, c1=(a1*a1+b1*b1)/2; 27 double a2=c.x-a.x, b2=c.y-a.y, c2=(a2*a2+b2*b2)/2; 28 double d=a1*b2-a2*b1; 29 struct node tmp; 30 tmp.x = a.x + (c1*b2-c2*b1)/d ; 31 tmp.y = a.y+(a1*c2-a2*c1)/d; 32 return tmp; 33 } 34 35 void min_cover_circle(){ 36 random_shuffle(p,p+n); 37 central = p[0]; 38 int i,j,k; 39 R=0; 40 for(i=1;i<n;i++){ 41 if(dist(central,p[i])+eps>R){ 42 central = p[i]; 43 R=0; 44 for(j=0;j<i;j++){ 45 if(dist(central,p[j])+eps>R){ 46 central.x=(p[i].x+p[j].x)/2; 47 central.y=(p[i].y+p[j].y)/2; 48 R=dist(central,p[j]); 49 for(k=0;k<j;k++){ 50 if(dist(central,p[k])+eps>R){ 51 central=circumcenter(p[i],p[j],p[k]); 52 R=dist(central,p[k]); 53 } 54 } 55 } 56 } 57 } 58 } 59 } 60 61 62 int main(){ 63 double x,y; 64 while(scanf("%lf%lf%d",&x,&y,&n)!=EOF){ 65 int i; 66 for(i=0;i<n;i++){ 67 scanf("%lf%lf",&p[i].x,&p[i].y); 68 } 69 min_cover_circle(); 70 printf("(%.1lf,%.1lf). %.1lf ",central.x,central.y,R); 71 } 72 return 0; 73 }
模拟退火算法
1 /* 2 参考：http://blog.sina.com.cn/s/blog_64675f540100sehz.html 3 题意：给定n个点，找到一个点，使得n个点到这个点的和值最小 4 模拟退火法 5 模拟退火的过程 6 1 找到这些点所在的范围，用两个点框定（代码中e1，e2两个点) 7 2 在这个范围内生成NUM个点（NUM自定） 8 3 对于每个生成的点i，在其周围生成NUM个点，一旦有点优于i，则替换。 9 4 缩小范围D，若D<精度，退出，否则执行 3 10 5 遍历所有NUM个点，找到val的最大值 11 */ 12 #include <iostream> 13 #include <cstdio> 14 #include <cstdlib> 15 #include <cmath> 16 #include <cstring> 17 #include <string> 18 #include <algorithm> 19 #include <set> 20 #include<map> 21 #include<ctime> 22 using namespace std; 23 const int NUM=30; 24 const int RAD=1000; 25 struct point 26 { 27 double x,y,val; 28 point(){} 29 point(double _x,double _y):x(_x),y(_y){} 30 }p[1001],May[NUM],e1,e2; 31 int n; 32 double dis(point a,point b) 33 { 34 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 35 } 36 double judge(point t)//评价函数，得到点t的评价值val 37 { 38 double len=0; 39 for(int i=0;i<n;i++) 40 len=max(len,dis(t,p[i])); 41 return len; 42 } 43 double Rand(){return rand()%(RAD+1)/(1.0*RAD);}//随机产生0-1的浮点数 44 point Rand_point(point a,point b)//在a,b框定的四边形内随机生成点 45 { 46 point tmp=point(a.x+(b.x-a.x)*Rand(),a.y+(b.y-a.y)*Rand()); 47 tmp.val=judge(tmp); 48 return tmp; 49 } 50 void solve(double D) 51 { 52 for(int i=0;i<NUM;i++) 53 May[i]=Rand_point(e1,e2);//步骤2 54 while(D>0.01)//步骤 3 55 { 56 for(int i=0;i<NUM;i++) 57 for(int j=0;j<NUM;j++) 58 { 59 point tmp=Rand_point(point(May[i].x-D,May[i].y-D),point(May[i].x+D,May[i].y+D)); 60 if(tmp.val<May[i].val) 61 { 62 May[i]=tmp; 63 } 64 } 65 D*=0.5; 66 } 67 point ans; 68 ans.val=1LL<<45; 69 for(int i=0;i<NUM;i++) 70 if(May[i].val<ans.val) 71 ans=May[i]; 72 printf("(%.1f,%.1f). ",ans.x,ans.y); 73 printf("%0.1f ",ans.val); 74 } 75 int main() 76 { 77 srand(time(0)); 78 e2=point(0,0); 79 while(scanf("%lf%lf%d",&e1.x,&e2.y,&n)!=EOF) 80 { 81 for(int i=0;i<n;i++) 82 { 83 scanf("%lf%lf",&p[i].x,&p[i].y); 84 e1.x=min(e1.x,p[i].x);//框定初始范围 85 e1.y=min(e1.y,p[i].y); 86 e2.x=max(e2.x,p[i].x); 87 e2.y=max(e2.y,p[i].y); 88 } 89 solve(max(e2.y-e1.y,e2.x-e1.x)); 90 } 91 }
2.最小矩形包围所有点（矩形的边平行于坐标轴）
1 #include <stdio.h> 2 #define ARRAY_SIZE 100 3 struct Point { 4 int x; 5 int y; 6 }; 7 struct Point p[ARRAY_SIZE]; 8 int main(void) { 9 int i, n; 10 int min_x, max_x, min_y, max_y; 11 scanf("%d", &n); 12 for(i = 0; i < n; ++i) 13 scanf("%d%d", &p[i].x, &p[i].y); 14 min_x = max_x = p[0].x; 15 min_y = max_y = p[0].y; 16 for(i = 1; i < n; ++i) { 17 if(min_x > p[i].x) 18 min_x = p[i].x; 19 if(max_x < p[i].x) 20 max_x = p[i].x; 21 if(min_y > p[i].y) 22 min_y = p[i].y; 23 if(max_y < p[i].y) 24 max_y = p[i].y; 25 } 26 printf("%d %d %d %d ", min_x, min_y, max_x, max_y); 27 return 0; 28 }
3.最小矩形包围所有点（矩形的边可以不平行于坐标轴）例题：2017_SWERC_K
1 //最小矩形包围点 2 #include <bits/stdc++.h> 3 typedef long long ll; 4 using namespace std; 5 //平面点的结构体模板 6 struct Point{ 7 double x,y; 8 Point(){} 9 Point(double _x,double _y){ x=_x; y=_y; } 10 Point operator + (Point p){ return Point(x+p.x,y+p.y); } 11 Point operator - (Point p){ return Point(x-p.x,y-p.y); } 12 Point operator * (double d){ return Point(d*x,d*y); } 13 double dot(Point p){ return x*p.x+y*p.y; } //内积 14 //外积，外积等于0则两点连线过原点或矢量平行，>0则连线斜率>45°，<0则连线斜率<45° 15 double det(Point p){ return x*p.y-p.x*y; } 16 }ps[200010]; 17 //加上操作变成凸包模板 18 int n; 19 bool cmpxy(const Point &p,const Point &q){ //排序 20 if(p.x!=q.x) return p.x<q.x; 21 return p.y<q.y; 22 } 23 vector<Point> convex_hull(Point *ps,int n){ 24 sort(ps,ps+n,cmpxy); 25 int k=0; 26 vector<Point> qs(n*2); //构造中的凸包 27 for(int i=0;i<n;++i){ //构造凸包的下侧 28 //相邻边的叉积大于0 29 while(k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0)--k; 30 qs[k++]=ps[i]; 31 } 32 for(int i=n-2,t=k;i>=0;--i){ //构造凸包的上侧 33 //相邻边的叉积小于0 34 while(k>t && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0)--k; 35 qs[k++]=ps[i]; 36 } 37 qs.resize(k-1); //去掉重复的左下角点 38 return qs; 39 } 40 double dist(Point p,Point q){ //两点距离平方 41 return (p-q).dot(p-q); 42 } 43 int main() 44 { 45 int n,r; 46 while(~scanf("%d%d",&n,&r)&&n) 47 { 48 for(int i=0;i<n;i++) 49 { 50 scanf("%lf%lf",&ps[i].x,&ps[i].y); 51 } 52 vector<Point> qs=convex_hull(ps,n); 53 int p=1; 54 double len; 55 int sz=qs.size(); 56 if(sz==2) 57 { 58 printf("0.000000000 "); 59 continue; 60 } 61 double minofmaxh=2.0*r; 62 qs[sz]=qs[0]; 63 for(int i=0;i<sz;i++) 64 { 65 len=sqrt(dist(qs[i],qs[i+1])); 66 while(true) 67 { 68 double h=fabs((qs[i]-qs[i+1]).det(qs[p]-qs[i+1]))/len; 69 if(p<sz-1) 70 { 71 double h1=fabs((qs[i]-qs[i+1]).det(qs[p+1]-qs[i+1]))/len; 72 if(h1>=h) 73 { 74 p++; 75 } 76 else 77 { 78 minofmaxh=min(minofmaxh,h); 79 break; 80 } 81 } 82 else 83 { 84 double h1=fabs((qs[i]-qs[i+1]).det(qs[0]-qs[i+1]))/len; 85 if(h1>=h) 86 { 87 p=0; 88 } 89 else 90 { 91 minofmaxh=min(minofmaxh,h); 92 break; 93 } 94 } 95 } 96 } 97 printf("%.10lf ",minofmaxh); 98 } 99 return 0; 100 }
4.ACM-ICPC 2018 沈阳赛区网络预赛 E. The cake is a lie
题意：给你n个点圆的圆心和相同的半径R，求至少覆盖s个小圆的最小大圆半径是多少
题解：这题的前身可以认为是POJ1981，POJ1981讲的是平面上有n个点单位圆最多能覆盖几个点，他是有O(n2logn)O(n2logn)的做法的，那我们回到这一题，我们先不考虑整圆，只考虑圆心，如果我们求出来的大圆能包含s个圆心那么我们让大圆的半径增加R，那就可以覆盖整个圆了，所以现在的问题就是如何快速的求出对一个大圆的半径有最多能覆盖多少个点。
附上POJ1981的详解链接
1 朴素算法：由限制条件可知，最后目标圆圆心一定可由某两点唯一确定（即目标圆的圆弧上一定恰有两个平面点集中的点）。
具体实现：C(N,2)枚举任意两圆，构造可行圆，枚举N个点，判断是否在圆上或者圆内（即是否被覆盖）。复杂度O(N^3)。我所知的标准做法：O(n^2*logn)的圆上交点扫描法。 2 转化为圆上的弧最多被覆盖几次的问题。 3 直觉可解，圆上某个弧如果被某个圆覆盖过，那么我们把目标圆圆心放在该段弧上时，目标圆一定能把该弧所在圆圆心覆盖（即平面点集中的某点）。 4 所以，被覆盖次数最多的弧的被覆盖次数即为目标圆（最优）最多能够覆盖的点数。 5 具体做法： 6 依旧是C(N,2)枚举任意两点，以BA为例，此时i=1,j=0。 7 求取BA的极角theta，利用atan2（因变量范围(-pi,pi]）。 8 为方便极角排序，对于负角要+2pi。 9 B、A构造的圆交于X1，X4两点，求取X1B与BA的夹角和X4B与AB的夹角（相同）。 10 Phi ，acos(D/2.0/R)，D=dis(A,B)。 11 以上可得出两圆交点在当前圆上对于圆心的极角。 12 为方便极角排序，每个极角均加2pi，使域扩张到[0，4pi]，从而防止了跨0角的出现。 13 alpha[t].value=theta – phi + 2pi,alpha[t].flag = 1(弧的进入端) 14 alpha[t+1].value = theta + phi + 2pi,alpha[t+1].flag = false（弧的退出端）。 15 接下来，对于每个圆上的t（偶数）个极角进行排序，再顺序扫描，动态更新最大值。 16 for(i,0,t) 17 if(alpha[i].flag) 18 sum++ 19 else 20 sum–— 21 if(sum > ans) 22 ans = sum 23 最后输出答案时，输出ans+1即可。（因为要算上弧所在圆自身）
模板如下：
#define Mn 300 const double eps = 1e-9; const double pi = acos(-1.0); #define sqr(x) ((x) * (x)) struct point { double x,y; void read() { scanf("%lf%lf",&x,&y); } void print() { printf("%lf%lf ",x,y); } double friend dis(const point &a,const point &b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)); } } p[Mn + 5]; struct alpha { double v; bool flag; bool friend operator <(const alpha &a,const alpha &b) { return a.v < b.v; } } alp[Mn * Mn + 5]; int solve(int n,double R)//半径为R的圆最多能覆盖多少个点 { int MAX = 0; for(int i = 0; i < n; i++) { int t = 0; for(int j = 0; j < n; j++) { if(i == j) continue; double theta,phi,D; D = dis(p[i],p[j]); if(D > 2.0 * R) continue; theta = atan2(p[j].y - p[i].y,p[j].x - p[i].x); if(theta < 0) theta += 2 * pi; phi = acos(D / (2.0 * R)); alp[t].v = theta - phi + 2 * pi; alp[t].flag = true; alp[t + 1].v = theta + phi + 2 * pi; alp[t + 1].flag = false; t += 2; } sort(alp,alp + t); int sum = 0; for(int j = 0; j < t; j++) { if(alp[j].flag) sum ++; else sum --; if(sum > MAX) MAX = sum; } } return MAX+1; }
题解来自dalao：https://blog.csdn.net/ftx456789/article/details/82532729（如有冒犯，请留言，我来删之QAQ，逃
本题网络赛E题解：二分答案，那么每次check相当于是给定一个半径的圆，然后问这个圆最多覆盖多少个点，我们可以枚举一个点，然后再枚举每个与他距离<=2r的点，就可以求出所有的相交弧，离散化之后，求出覆盖最多次的弧，就是答案了。复杂度O(n2log(n)⋅log(30000))O(n2log(n)⋅log(30000))。
1 #include <bits/stdc++.h> 2 using namespace std; 3 const int maxn = 605; 4 const double eps = 1e-6; 5 int i, j, t, n, m, k, z, y, x; 6 double l, r, mid; 7 int R; 8 int cmp(double x) 9 { 10 if (fabs(x) < eps) return 0; 11 if (x > 0) return 1; 12 return -1; 13 } 14 struct point 15 { 16 double x, y; 17 point() {} 18 point(double _x, double _y) 19 : x(_x), y(_y) {} 20 void input() 21 { 22 scanf("%lf%lf", &x, &y); 23 } 24 friend point operator+(const point &a, const point &b) 25 { 26 return point(a.x + b.x, a.y + b.y); 27 } 28 friend point operator-(const point &a, const point &b) 29 { 30 return point(a.x - b.x, a.y - b.y); 31 } 32 double norm() 33 { 34 return sqrt(x * x + y * y); 35 } 36 double angl() 37 { 38 return atan2(y, x); 39 } 40 } poi[maxn]; 41 double dis[maxn][maxn], ang[maxn][maxn]; 42 pair<double, int> pdi[maxn]; 43 #define mp(a, b) make_pair(a, b) 44 bool check(double r) 45 { 46 int i, j, t, ans = 1, k; 47 double d; 48 for (i = 1; i <= n; i++) 49 { 50 t = 0; 51 for (j = 1; j <= n; j++) 52 if (i != j) 53 { 54 if (cmp(dis[i][j] - 2.0 * r) > 0) continue; 55 d = acos(dis[i][j] / (2.0 * r)); 56 pdi[++t] = mp(ang[i][j] - d, 1); 57 pdi[++t] = mp(ang[i][j] + d, -1); 58 } 59 sort(pdi + 1, pdi + t + 1); 60 k = 1; 61 for (j = 1; j <= t; j++) 62 { 63 k += pdi[j].second; 64 ans = max(ans, k); 65 } 66 } 67 return ans >= m; 68 } 69 int main() 70 { 71 // cout << time(NULL) << endl; 72 // freopen("tcial.in", "r", stdin); 73 // freopen("tcial.out", "w", stdout); 74 int T, I; 75 scanf("%d", &T); 76 for (I = 1; I <= T; I++) 77 { 78 scanf("%d%d", &n, &m); 79 for (i = 1; i <= n; i++) poi[i].input(); 80 scanf("%d", &R); 81 if (m > n) 82 { 83 printf("The cake is a lie. "); 84 continue; 85 } 86 for (i = 1; i <= n; i++) 87 for (j = 1; j <= n; j++) 88 dis[i][j] = (poi[j] - poi[i]).norm(), ang[i][j] = (poi[j] - poi[i]).angl(); 89 l = R; 90 r = 30000; 91 while (r - l > eps) 92 { 93 mid = (l + r) / 2; 94 if (check(mid - R)) 95 r = mid; 96 else 97 l = mid; 98 } 99 printf("%.6f ", r); 100 } 101 // cout << time(NULL) << endl; 102 return 0; 103 }
5.最小覆盖圆衍生问题：如何求能覆盖N个点中M个点(N≥M)的最小圆？
详情请见知乎上dalao的畅谈，菜鸡膜拜。链接