\(\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);

}

巴特西

JZOJ 1082. 【GDOI2005】选址

\(\text{Problem}\)

\(\text{Solution}\)

\(\text{Code}\)

最新文章

热门文章