JZOJ 1082. 【GDOI2005】选址

\(\text{Problem}\)

很久以前，在世界的某处有一个形状为凸多边形的小岛，岛上的居民们决定建一个祭坛，居民们任务祭坛的位置离岛的顶点处越远越好。
你的任务是求凸多边形内一点，使其与各顶点的距离中最短的距离最远，点在边上也可以。 这样的点可能有多个，你只需输出这些点与各顶点的最短距离。

\(\text{Solution}\)

非常经典的题
以答案为半径做圆，满足圆心到所有点距离大于等于半径
考虑二分半径，判断圆心存不存在
枚举任意两点，考虑分别以此半径做圆交出的点（选取内部的点，叉积判断是否在内部）
如果这点没有被所有圆覆盖，即这点到凸边形顶点的最短距离大于等于半径，说明圆心存在
算交点用相似，如果没有交点考虑其在边上交出的点一样判断
精度很神奇，不要用 \(\text{long double}\)
本地过不了数据1 \(OJ\) 上却过了？！

\(\text{Code}\)

#include <cstdio>
#include <algorithm>
#include <cmath>
#define RE register
#define IN inline
using namespace std;

const int N = 105;
const double eps = 1e-8;
double area;
int n;
struct Vector{
	double x, y;
	IN Vector(double xx = 0, double yy = 0){x = xx, y = yy;}
	IN Vector operator - (const Vector &B){return Vector(x - B.x, y - B.y);}
	IN double operator * (const Vector &B){return fabs(x * B.y - y * B.x);}
}p[N];

IN double sqr(double x){return x * x;}
IN double distance(Vector a, Vector b){return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));}
IN double Area(Vector a, Vector b, Vector c){return (a - b) * (c - b);}
IN double sum_Area(Vector a)
{
	double res = 0;
	for(RE int i = 1; i < n; i++) res += Area(a, p[i], p[i + 1]);
	return res + Area(a, p[n], p[1]);
}
IN int isIn(Vector a){if (fabs(sum_Area(a) - area) <= eps) return 1; return 0;}
IN int OK(Vector a, double r)
{
	double res = 1e18;
	for(RE int i = 1; i <= n; i++) res = min(res, distance(a, p[i]));
	return res >= r;
}

IN int check(double r)
{
	for(RE int i = 1; i <= n; i++)
		for(RE int j = 1; j < i; j++)
		{
			double d = distance(p[i], p[j]);
			if (d > r * 2)
			{
				double k = (p[i].y - p[j].y) / (p[i].x - p[j].x), x, y;
				y = min(p[i].y, p[j].y) + fabs(p[i].y - p[j].y) * r / d;
				if (k > 0)
				{
					x = min(p[i].x, p[j].x) + fabs(p[i].x - p[j].x) * r / d;
					if (OK(Vector{x, y}, r)) return 1;
					x = p[i].x + p[j].x - x, y = p[i].y + p[j].y - y;
					if (OK(Vector{x, y}, r)) return 1;
				}
				else{
					x = max(p[i].x, p[j].x) - fabs(p[i].x - p[j].x) * r / d;
					if (OK(Vector{x, y}, r)) return 1;
					x = p[i].x + p[j].x - x, y = p[i].y + p[j].y - y;
					if (OK(Vector{x, y}, r)) return 1;
				}
			}
			else{
				double w = sqrt(r * r - d * d / 4), k = w / d;
				double dx = (p[i].x + p[j].x) / 2, dy = (p[i].y + p[j].y) / 2;
				double x = dx - fabs(p[i].y - p[j].y) * k, y = dy + fabs(p[i].x - p[j].x) * k;
				if (isIn(Vector{x, y}) && OK(Vector{x, y}, r)) return 1;
				x = p[i].x + p[j].x - x, y = p[i].y + p[j].y - y;
				if (isIn(Vector{x, y}) && OK(Vector{x, y}, r)) return 1;
			}
		}
	return 0;
}

int main()
{
	scanf("%d", &n);
	for(RE int i = 1; i <= n; i++) scanf("%lf%lf", &p[i].x, &p[i].y);
	area = sum_Area(p[1]); double r = 0;
	for(RE int i = 1; i <= n; i++)
		for(RE int j = 1; j < i; j++) r = max(r, distance(p[i], p[j]));
	double l = 0, mid = (l + r) / 2, ans;
	for(RE int i = 0; i < 60; i++, mid = (l + r) / 2)
		if (check(mid)) ans = mid, l = mid; else r = mid;
	printf("%.3lf\n", ans);
}