题解 P2533 【[AHOI2012]信号塔】

· · 题解

原题

上次的题解被学长hack掉了，懒得调了，把退火做法放上来吧。

考虑模拟退火：选一个初始点，从这个点开始跑退火。问题是初始点怎么选。

首先所有点取平均大概是不怎么对的，只要很多点密集分布在一个区域再放一个点在很远处就可以卡掉。

发现这个圆只要能够盖住这些点的凸包就能够盖住所有点，所以考虑在凸包内找一个比较中间的点。比方说凸包直径的中点。

于是把凸包直径的中点用旋转卡壳卡出来，用该点作为初始点，退火随机偏移，记录最优解。

当然你也可以选择其他较好的点作为初始点，这样可以省去旋转卡壳的代码。然而我没有写出来。

代码

const int N = 1e6 + 10;
const int inf = 0x3fffffff;
const double eps = 1e-6;
int n;
struct node
{
    double x,y;
    friend bool operator == (node x , node y)
        {return x.x == y.x && x.y == y.y;}
}a[N],*st;
int s[N],top;
double ct(node p1 , node p2) // 叉积 
{
    p1.x -= st->x,p2.x -= st->x;
    p1.y -= st->y,p2.y -= st->y;
    return p1.x * p2.y - p2.x * p1.y;
}
inline double dis(node p1 , node p2) // 距离 
    {return sqrt( (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) );}
inline bool cmp(node p1 , node p2) // 极角排序 
    {return atan2(p1.y - st->y , p1.x - st->x) < atan2(p2.y - st->y , p2.x - st->x);}
inline node mid(node p1 , node p2)
    {return (node){ (p1.x + p2.x) / 2 , (p1.y + p2.y) / 2 };}
inline void graham() // Graham凸包 
{
    st = &a[1];
    for(register int i = 1;i <= n;i++)
        if(a[i].y < st->y || a[i].y == st->y && a[i].x < st->x)
            st = &a[i];
    swap(a[1] , *st);
    st = &a[1];
    sort(a + 2 , a + n + 1 , cmp);
    top = 1,s[0] = 1,s[1] = 2;
    for(int i = 3;i <= n;i++)
    {
        st = &a[ s[top - 1] ];
        while(top > 0 && ct( a[ s[top] ] , a[i] ) <= 0)
            st = &a[ s[--top - 1] ];
        s[++top] = i;
    }
}
double query(double x , double y) // 以(x, y)为圆心时的最小半径 
{
    double res = -inf;
    for(register int i = 0;i <= top - 1;i++)
        res = max(dis( (node){x , y} , a[ s[i] ] ) , res);
    return res;
}
const double delta = 0.997; // 降温系数 
double ans,ansx,ansy;
void SA() // 模拟退火 
{
    double T = 3000;
    while(T > 1e-12)
    {
        double x = ansx + (rand() * 2 - RAND_MAX) * T;
        double y = ansy + (rand() * 2 - RAND_MAX) * T;
        double cur = query(x , y);
        if(cur < ans)
            ans = cur,ansx = x,ansy = y;
        else if( exp( (cur - ans) / T ) * RAND_MAX < rand() ) // 一定概率接受当前答案 
            ansx = x,ansy = y;
        T *= delta;
    }
}
int main()
{
    srand( time(NULL) );
    scanf("%d" , &n);
    for(int i = 1;i <= n;i++)
        scanf("%lf%lf" , &a[i].x , &a[i].y);
    graham();
    s[++top] = s[0];
    node A,B; // 凸包直径上的两点 
    int j = 1;
    double maxn = -inf;
    for(register int i = 0;i <= top - 1;i++) // 旋转卡壳
    {
        st = &a[ s[i] ];
        while(fabs( ct( a[ s[i + 1] ] , a[ s[j] ] ) ) - fabs( ct( a[ s[i + 1] ] , a[ s[j + 1] ] ) ) < eps)
        {
            ++j;
            if(j >= top)
                j = 0;
        }
        double d = dis( a[ s[i] ] , a[ s[j] ] );
        if(d > maxn)
        {
            maxn = d;
            A = a[ s[i] ];
            B = a[ s[j] ];
        }
    }
    ansx = (A.x + B.x) / 2;
    ansy = (A.y + B.y) / 2;
    ans = query(ansx , ansy);
//  while(clock() < CLOCKS_PER_SEC * 0.95)
    for(int i = 1;i <= 150;i++)
        SA();
    printf("%.2lf %.2lf %.2lf\n" , ansx , ansy , ans);
    return 0;
}