随机增量法：bzoj 1336 bzoj 1337 最小圆覆盖

1337: 最小圆覆盖

Time Limit: 1 Sec Memory Limit: 64 MB
Submit: 1170 Solved: 573
[Submit][Status][Discuss]

Description

给出平面上N个点,N<=10^5.请求出一个半径最小的圆覆盖住所有的点

Input

第一行给出数字N，现在N行，每行两个实数x,y表示其坐标.

Output

输出最小半径，输出保留三位小数.

Sample Input

1 0

0 1

0 -1

-1 0

Sample Output

1.000

题目：求出n个点的最小圆覆盖

随机增量法：

①将数组重新打乱（random_shuffle()函数）

②假设已经求出了前i个点的最小圆覆盖为Ci，检查第i+1个点是否在圆内，如果在，则C(i+1)=Ci并继续判定下一个点，不在的话，就令C(i+1)的圆心为第1个点和第i+1个点连线的中点，半径为连线长度的一半，并执行步骤③

③重新检查前i个点是否在圆C(i+1)上，如果某个点j不在，就用类似步骤②的方法继续修改圆心和半径使得j点在圆C(i+1)上，并执行步骤④

④再来一次一模一样的检查过程，不过这次只要有点k不在圆上，就令圆C(i+1)的圆心为三角形(i, j, k)的外心并扩增圆C(i+1)的半径

搞定，复杂度O(n^3)不过平均复杂度只有O(n)，证明略，应该很容易找到证明

#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
#define eps 1e-8
typedef struct Point
{double x;double y;
}Point;
Point s[100005];
Point operator + (const Point a, const Point b)
{  Point k;k.x = a.x+b.x, k.y = a.y+b.y;return k;
}
Point operator - (const Point a, const Point b)
{Point k;k.x = a.x-b.x, k.y = a.y-b.y;return k;
}
Point operator / (const Point a, double x)
{Point k;k.x = a.x/x, k.y = a.y/x;return k;
}
Point operator * (const Point a, double x)
{Point k;k.x = a.x*x, k.y = a.y*x;return k;
}double Dist(Point a, Point b)
{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point Waixin(Point a, Point b, Point c)     //求三角形外接圆圆心（外心）
{Point ans;double a1, a2, b1, b2, c1, c2;a1 = 2*(b.x-a.x), b1 = 2*(b.y-a.y); c1 = (b.x*b.x)-(a.x*a.x)+(b.y*b.y)-(a.y*a.y);a2 = 2*(c.x-a.x), b2 = 2*(c.y-a.y); c2 = (c.x*c.x)-(a.x*a.x)+(c.y*c.y)-(a.y*a.y);if(fabs(a1)<eps)ans.y = c1/b1, ans.x = (c2-ans.y*b2)/a2;else if(fabs(b1)<eps)ans.x = c1/a1, ans.y = (c2-ans.x*a2)/b2;elseans.x = (c2*b1-c1*b2)/(a2*b1-a1*b2), ans.y = (c2*a1-c1*a2)/(b2*a1-b1*a2);return ans;
}
int main(void)
{Point P;double R;int n, i, j, k;scanf("%d", &n);for(i=1;i<=n;i++)scanf("%lf%lf", &s[i].x, &s[i].y);random_shuffle(s+1, s+n+1);P = s[1], R = 0;for(i=2;i<=n;i++){if(Dist(P, s[i])-R>eps){P = (s[1]+s[i])/2;R = Dist(P, s[i]);for(j=2;j<=i-1;j++){if(Dist(P, s[j])-R>eps){P = (s[i]+s[j])/2;R = Dist(P, s[i]);for(k=1;k<=j-1;k++){if(Dist(P, s[k])-R>eps){P = Waixin(s[i], s[j], s[k]);R = Dist(P, s[i]);}}}}}}printf("%.3f\n", R);return 0;
}