solution-p2533

题解 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;
}