bzoj1337 minimum circle cover

Description

Give N points on the plane, N < = 10 ^ 5. Request a circle with the smallest radius to cover all points

Solution

Of course, new things need porcelain. I can't remember the computational geometry very much. What's more, it's such a magical algorithm
Obviously, an optimal circle must pass at least two points, so we default 1 as the initial center.
Every time the i-th point is added, if not all points can be covered, it means that the i-th point must be on the circle. Go back to enumerate another point to make a new circle and judge again
If the circle made by two points still doesn't satisfy the full coverage, then find the third point so that the three points determine a circle. This circle is the answer of the first i points

In order to ensure the complexity, n points need to be randomly arranged. Since the probability outside the circle is about 3/i by adding a point to i points, the expectation of time complexity is O (n)

I'm probably a mentally retarded player. Push the content of middle school for half a day to find the intersection point of linear equation, gg

Code

#include <cstdio>
#include <cmath>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
const int N=200005;
struct pos {double x,y;}p[N],cir;
double r=0;
double sqr(double x) {return x*x;}
double get_dis(pos a,pos b) {return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}
pos get_mid(pos a,pos b) {return (pos){(a.x+b.x)*0.5,(a.y+b.y)*0.5};}
pos get_center(pos a,pos b,pos c) {
    pos cen1=get_mid(a,b);
    pos cen2=get_mid(a,c);
    pos ret;
    double k1,k2,b1,b2;
    if (a.y==b.y) {
        k2=-1.0/((c.y-a.y)/(c.x-a.x));
        b2=cen2.y-k2*cen2.x;
        ret=(pos){cen1.x,cen1.x*k2+b2};
    } else if (a.y==c.y) {
        k1=-1.0/((b.y-a.y)/(b.x-a.x));
        b1=cen1.y-k1*cen1.x;
        ret=(pos){cen2.x,cen2.x*k1+b1};
    } else {
        k1=-1.0/((b.y-a.y)/(b.x-a.x));
        b1=cen1.y-k1*cen1.x;
        k2=-1.0/((c.y-a.y)/(c.x-a.x));
        b2=cen2.y-k2*cen2.x;
        ret.x=(b1-b2)/(k2-k1); ret.y=k1*ret.x+b1;
    }
    return ret;
}
int main(void) {
    int n; scanf("%d",&n);
    rep(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y);
    std:: random_shuffle(p+1,p+n+1);
    cir=p[1];
    rep(i,2,n) {
        double d=get_dis(p[i],cir);
        if (d<=r) continue;
        r=0; cir=p[i];
        rep(j,1,i-1) {
            if (get_dis(p[j],cir)<=r)  continue;
            r=get_dis(p[j],p[i])*0.5;
            cir=get_mid(p[j],p[i]);
            rep(k,1,j-1) {
                if (get_dis(p[k],cir)<=r) continue;
                cir=get_center(p[i],p[j],p[k]);
                r=get_dis(cir,p[i]);
            }
        }
    }
    printf("%.3lf\n", r);
    return 0;
}

Posted on Mon, 04 May 2020 11:43:07 -0700 by kubis