[bzoj4138][随机化]最小覆盖双圆问题

4138: [FJOI2015]最小覆盖双圆问题

Time Limit: 10 Sec Memory Limit: 128 MB
Submit: 171 Solved: 42
[Submit][Status][Discuss]
Description

给定平面上n个点(x1,y1),…,(xn,yn),找出2个半径相同的圆R1和R2,覆盖给定的n个点,且半径最小。

设计一个算法,计算出所求最小覆盖双圆 R1 和 R2 的半径。
Input

输入有多个测试实例。每个实例的第1行中给出正整数n,n<1000,表示平面上有n个点。
接下来的n行中每行给出2个实数(x, y),-100000≤x≤100000,-100000≤y≤100000。
最后一行有一个0表示结束。
Output

对于每组数据,输出最小的符合题意的圆的半径,保留两位小数。
Sample Input

3

0.00 0.00

1.00 0.00

0.00 4.00

10

0.00 0.00

0.00 3.00

1.00 6.00

2.00 2.00

3.00 5.00

5.00 3.00

6.00 3.00

9.00 5.00

10.00 5.00

11.00 3.00

0
Sample Output

0.50

3.05

HINT

对于100%的数据,n<=1000,|xi|,|yi|<=100000,(T<=10)

Source

sol:

人类智慧题
做法T了,大概60分
考虑有一条直线和y平行,用这条直线来划分点集,那么这条直线右移的时候肯定能使得R1变大,R2变小,当R1>R2后这条直线就到了一个合法的位置了。2个点集就是各自做单圆覆盖,ans=min(max)。但是感觉就会出问题对吧。我们把坐标系随机转一下,发现可以过很多点了。。

#include<iostream>
#include<algorithm>
#include<string>
#include<cstring>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<ctime>
using namespace std;
typedef long long ll;
typedef double s64;
int n,m;
const int N=1100;
const s64 pi=acos(-1);
const s64 eps=1e-7;
s64 cc,dd;
struct P
{
    s64 x,y;
    friend P operator -(P a,P b) {return (P){a.x-b.x,a.y-b.y};}
    friend P operator +(P a,P b) {return (P){b.x+a.x,b.y+a.y};}
    friend s64 operator *(P a,P b) {return a.x*b.y-a.y*b.x;}
    friend P operator *(P a,s64 b) {return (P){a.x*b,a.y*b};}
    friend s64 operator /(P a,P b) {return b.x*a.x+b.y*a.y;}
    friend bool operator <(P a,P b) {return a.x<b.x||a.x==b.x&&a.y<b.y;}
    friend P rev(P a) {return (P){-a.y,a.x};}
    friend P rotate(P &a,s64 b)
    {
        return (P){a.x*cc-a.y*dd,a.x*dd+a.y*cc};
    }
}a[N];
struct L
{
    P a,b,v;
    friend inline P inter(const L &a,const L &b)
    {
        P nw=b.a-a.a;
        s64 tt=(nw*a.v)/(a.v*b.v);
        return b.a+b.v*tt;
    }
};
inline s64 sqr(const s64 &x)
{
    return x*x;
}
inline s64 dis(const P &a,const P &b)
{
    return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
inline P find(P a,P b,P c)
{
    if(fabs((b-a)*(c-b))<=eps)
    {
        if(b.x>c.x) swap(b,c);
        if(a.x>b.x) swap(a,b);
        if(b.x>c.x) swap(b,c);
        return (a+c)*0.5;
    }
    L x,y;
    x.a=(a+b)*0.5;
    x.v=rev(b-a);
    x.b=x.a+x.v;
    y.a=(b+c)*0.5;
    y.v=rev(c-b);
    y.b=y.a+y.v;
    return inter(x,y);
}
inline s64 cover(P *a,int n)
{
    random_shuffle(a+1,a+1+n);
    P t; s64 r;
    t=a[1];r=0.0;
    for(int i=2;i<=n;++i)
    if(dis(t,a[i])>r+eps)
    {
        t=(a[1]+a[i])*0.5;
        r=dis(a[1],t);
        for(int j=2;j<i;++j)
        if(dis(t,a[j])>r+eps)
        {
            t=(a[i]+a[j])*0.5;
            r=dis(a[i],t);
            for(int k=1;k<j;++k)
            if(dis(t,a[k])>r+eps)
            {
                t=find(a[i],a[j],a[k]);
                r=dis(t,a[i]);
            }
        }
    }
    return r;
}
P b[N],c[N],d[N];
s64 ans;
inline bool check(int n)
{
    for(int i=1;i<=n;++i) b[i]=a[i];
    for(int i=n+1;i<=m;++i) c[i-n]=a[i];
    s64 fx,fy;
    fx=cover(b,n);
    fy=cover(c,m-n);
    ans=min(ans,max(fx,fy));
    if(fx<=fy) return 1;
    return 0;
}
int main()
{
 //   freopen("data.in","r",stdin);
 //   freopen("data.out","w",stdout);
    srand(3268);
    while(true)
    {
        ans=1e9;
        cin>>m;
        if(!m) break;
        for(int i=1;i<=m;++i)
        scanf("%lf%lf",&d[i].x,&d[i].y);
        for(int o=1;o<=200;++o)
        {
            s64 th=(s64)rand()/32768.0*pi;
            cc=cos(th);
            dd=sin(th);
            for(int i=1;i<=m;++i)
            a[i]=rotate(d[i],th);
            sort(a+1,a+1+m);
            int l=1,r=m,mid;
            while(l+1<r)
            {
                mid=l+r>>1;
                if(check(mid)) l=mid;
                else r=mid;
            }
            check(l);
            if(l+1<=m) check(l+1);
            check(r);
        }
        printf("%.2lf\n",ans);
    }
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值