题意是，给出一些直径为1的树的位置，要求求出最小的能够把所有的树围起来的圆的半径。
开始的时候打算用类似模拟退火的方法，假定一个位置是最小圆的圆心，然后向每一个方向尝试移动圆心。如果能得到更小的一个圆，那么就将圆心移动到这里。如果把所有的方向都找遍都不能找到更小的圆，这时就要将点移动的距离减少，继续移动。直到移动的距离足够小的时候停止搜索。不过这样做，精度不能保证，所以最后还是要换回最小包围圈算法来通过这题。
下面的是包围圈暴力算法的代码，复杂度O(n^3)：
1 #include2 #include 3 #include 4 #include 5 #include 6 7 using namespace std; 8 9 const double PI = acos(-1.0); 10 const double EPS = 1e-8; 11 inline int sgn(double x) { return (EPS < x) - (x < -EPS);} 12 struct Point { 13 double x, y; 14 Point() {} 15 Point(double x, double y) : x(x), y(y) {} 16 bool operator < (Point a) const { return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;} 17 bool operator == (Point a) const { return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;} 18 } ; 19 typedef Point Vec; 20 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);} 21 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);} 22 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);} 23 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);} 24 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;} 25 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;} 26 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);} 27 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));} 28 inline double ptDis(Point a, Point b) { return vecLen(a - b);} 29 inline Vec normal(Vec x) { 30 double len = vecLen(x); 31 return Vec(-x.y / len, x.x / len); 32 } 33 inline Vec vecUnit(Vec x) { return x / vecLen(x);} 34 35 struct Line { 36 Point s, t; 37 Line() {} 38 Line(Point s, Point t) : s(s), t(t) {} 39 Point point(double x) { return s + (t - s) * x;} 40 } ; 41 typedef Line Seg; 42 Point lineIntersect(Point P, Vec v, Point Q, Vec w) { 43 Vec u = P - Q; 44 double t = crossDet(w, u) / crossDet(v, w); 45 return P + v * t; 46 } 47 inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);} 48 49 Line midLine(Point a, Point b) { 50 Point mid = (a + b) / 2.0; 51 return Line(mid, mid + normal(a - b)); 52 } 53 54 inline Point getPoint(Point a, Point b, Point c) { return lineIntersect(midLine(a, b), midLine(b, c));} 55 56 double getDis(Point a, Point b, Point c) { 57 if (sgn(crossDet(a, b, c)) == 0) return 0.0; 58 if (sgn(dotDet(a - b, c - b)) <= 0) return ptDis(a, c) / 2.0; 59 if (sgn(dotDet(a - c, b - c)) <= 0) return ptDis(a, b) / 2.0; 60 if (sgn(dotDet(b - a, c - a)) <= 0) return ptDis(b, c) / 2.0; 61 Point cen = getPoint(a, b, c); 62 if (sgn(ptDis(cen, a) - ptDis(cen, b)) || sgn(ptDis(cen, b) - ptDis(cen, c))) { 63 cout << "shit!" << endl; 64 while (1) ; 65 } 66 return ptDis(cen, a); 67 } 68 69 int andrew(Point *pt, int n, Point *ch) { 70 sort(pt, pt + n); 71 int m = 0; 72 for (int i = 0; i < n; i++) { 73 while (m > 1 && crossDet(ch[m - 2], ch[m - 1], pt[i]) <= 0) m--; 74 ch[m++] = pt[i]; 75 } 76 int k = m; 77 for (int i = n - 2; i >= 0; i--) { 78 while (m > k && crossDet(ch[m - 2], ch[m - 1], pt[i]) <= 0) m--; 79 ch[m++] = pt[i]; 80 } 81 if (n > 1) m--; 82 return m; 83 } 84 85 const int N = 111; 86 87 int main() { 88 // freopen("in", "r", stdin); 89 int n; 90 Point pt[N], ch[N]; 91 while (cin >> n && n) { 92 for (int i = 0; i < n; i++) cin >> pt[i].x >> pt[i].y; 93 n = andrew(pt, n, ch); 94 // cout << n << endl; 95 // for (int i = 0; i < n; i++) cout << ch[i].x << ' ' << ch[i].y << endl; 96 double dis = 0.0; 97 if (n == 2) { 98 dis = ptDis(ch[0], ch[1]) / 2.0; 99 } else {100 for (int i = 0; i < n; i++) {101 for (int j = i + 1; j < n; j++) {102 for (int k = j + 1; k < n; k++) {103 dis = max(dis, getDis(ch[i], ch[j], ch[k]));104 }105 }106 }107 }108 printf("%.2f\n", dis + 0.5);109 }110 return 0;111 }
之后将更新增量法的最小包围圈算法，平均复杂度O(n)。
增量法：
1 #include2 #include 3 #include 4 #include 5 #include 6 7 using namespace std; 8 9 const double EPS = 1e-10; 10 inline int sgn(double x) { return (x > EPS) - (x < -EPS);} 11 struct Point { 12 double x, y; 13 Point() {} 14 Point(double x, double y) : x(x),y(y) {} 15 bool operator < (Point a) const { return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;} 16 bool operator == (Point a) const { return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;} 17 Point operator + (Point a) const { return Point(x + a.x, y + a.y);} 18 Point operator - (Point a) const { return Point(x - a.x, y - a.y);} 19 Point operator * (double p) const { return Point(x * p, y * p);} 20 Point operator / (double p) const { return Point(x / p, y / p);} 21 } ; 22 typedef Point Vec; 23 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;} 24 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);} 25 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;} 26 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));} 27 inline Point normal(Vec x) { return Point(-x.y, x.x) / vecLen(x);} 28 29 Point lineIntersect(Point P, Vec v, Point Q, Vec w) { 30 Vec u = P - Q; 31 double t = crossDet(w, u) / crossDet(v, w); 32 return P + v * t; 33 } 34 inline Point getMid(Point a, Point b) { return (a + b) / 2.0;} 35 36 struct Circle { 37 Point c; 38 double r; 39 Circle() {} 40 Circle(Point c, double r) : c(c), r(r) {} 41 } ; 42 43 Circle getCircle(Point a, Point b, Point c) { 44 Vec v1 = b - a, v2 = c - a; 45 if (sgn(dotDet(b - a, c - a)) <= 0) return Circle(getMid(b, c), vecLen(b - c) / 2.0); 46 if (sgn(dotDet(a - b, c - b)) <= 0) return Circle(getMid(a, c), vecLen(a - c) / 2.0); 47 if (sgn(dotDet(a - c, b - c)) <= 0) return Circle(getMid(a, b), vecLen(a - b) / 2.0); 48 Point ip = lineIntersect(getMid(a, b), normal(v1), getMid(a, c), normal(v2)); 49 return Circle(ip, vecLen(ip - a)); 50 } 51 52 int andrew(Point *pt, int n, Point *ch) { 53 sort(pt, pt + n); 54 int m = 0; 55 for (int i = 0; i < n; i++) { 56 while (m > 1 && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--; 57 ch[m++] = pt[i]; 58 } 59 int k = m; 60 for (int i = n - 2; i >= 0; i--) { 61 while (m > k && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--; 62 ch[m++] = pt[i]; 63 } 64 if (n > 1) m--; 65 return m; 66 } 67 68 const int N = 555; 69 Point pt[N], ch[N]; 70 int rnd[N]; 71 72 void randPoint(Point *pt, int n) { 73 for (int i = 0; i < n; i++) rnd[i] = (rand() % n + n) % n; 74 for (int i = 0; i < n; i++) swap(pt[i], pt[rnd[i]]); 75 } 76 77 inline bool inCircle(Point p, Circle C) { return sgn(vecLen(C.c - p) - C.r) <= 0;} 78 79 int main() { 80 // freopen("in", "r", stdin); 81 int n; 82 while (cin >> n && n) { 83 for (int i = 0; i < n; i++) scanf("%lf%lf", &pt[i].x, &pt[i].y); 84 n = andrew(pt, n, ch); 85 randPoint(ch, n); 86 Circle ans = Circle(ch[0], 0.0), tmp; 87 for (int i = 0; i < n; i++) { 88 if (inCircle(ch[i], ans)) continue; 89 ans = Circle(ch[i], 0.0); 90 for (int j = 0; j < i; j++) { 91 if (inCircle(ch[j], ans)) continue; 92 ans = Circle(getMid(ch[i], ch[j]), vecLen(ch[i] - ch[j]) / 2.0); 93 for (int k = 0; k < j; k++) { 94 if (inCircle(ch[k], ans)) continue; 95 ans = getCircle(ch[i], ch[j], ch[k]); 96 } 97 } 98 } 99 // printf("%.2f %.2f %.2f\n", ans.c.x, ans.c.y, ans.r);100 printf("%.2f\n", ans.r + 0.5);101 }102 return 0;103 }
——written by Lyon