[LOJ#6437][BZOJ5373]「PKUSC2018」PKUSC

[LOJ#6437][BZOJ5373]「PKUSC2018」PKUSC

试题描述

九条可怜是一个爱玩游戏的女孩子。

最近她在玩一个无双割草类的游戏，平面上有 \(n\) 个敌人，每一个敌人的坐标为 \(x_i,y_i\)。可怜有一个技能是在平面上画一个 \(m\) 个点的简单多边形，并消灭所有严格在多边形内部的敌人。

不难发现如果想要快速的消灭敌人的话，只要画一个足够大的简单多边形就行了。但是这样的游戏性就太差了。于是可怜打算为游戏增加一定的随机性。

可怜在平面上随便画了一个 \(m\) 个点的简单多边形 \((a_i,b_i)\)。接下来可怜打算按照 \([-\pi,\pi)\) 上的均匀分布随机选取数字 \(\alpha\)（可以理解为等概率选取），并把这个简单多边形绕原点逆时针旋转 \(\alpha\) 的角度（弧度制）。

现在可怜给你了每一个点的坐标，多边形的坐标，你的任务是帮助可怜计算在随机旋转后她期望可以消灭多少个敌人。

输入

第一行四个整数 \(n,m\)。

接下来 \(n\) 行每行两个整数 \(x_i,y_i\) 描述了一个敌人的坐标。

接下来 \(m\) 行每行两个整数 \(a_i,b_i\) 按照逆时针的顺序描述了简单多边形的每一个顶点。

输出

输出一行一个整数表示期望值，保留五位小数。同时保证所有数据的小数点后第 \(6\) 位在舍入前不会是 \(4\) 和 \(5\)。

输入示例

4 4
0 0
1 0
-1 -1
0 1
0 0
1 0
1 1
0 1

输出示例

0.50000

数据规模及约定

对于 \(30\%\) 的数据，\(n,m \leq 15\)。

对于另外 \(30\%\) 的数据，选择的简单多边形是一个凸多边形。

对于 \(100\%\) 的数据，\(n \leq 200, m \leq 500, |x|,|y| \leq 10^6\).

题解

一道很裸的计算几何题，写一写练习一下，毕竟考场上这题爆零了很不爽……

考虑期望的线性性。我们可以对于每个敌人计算他被消灭的概率，将所有的概率累加即可。

多边形转可以变成固定多边形，然后转敌人的坐标。这样就变成了 \(n\) 次求圆有多少部分在多边形内部的问题。暴力求出圆与多边形的交点，然后对于每两个交点之间的弧，判一下这个弧的中点是否在多边形内部，如果在内部就将这个弧对应的圆心角 \(\theta\) 除以 \(2 \pi\) 的值累加即可。

需要的操作有：线段和圆求交点、线段和射线求交点（射线法）。

这个题需要调调 eps。eps 太小会 WA……

#include <bits/stdc++.h>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 510
#define nxt(x) ((x + 1) % n)
#define pre(x) ((x + n - 1) % n)

const double eps = 1e-7, pi = acos(-1.0);

int dcmp(double x) { return x < -eps ? -1 : (x > eps ? 1 : 0); }

struct Vec {
	double x, y;
	Vec(double _ = 0, double __ = 0): x(_), y(__) {}
	Vec operator + (const Vec &t) const { return Vec(x + t.x, y + t.y); }
	Vec operator - (const Vec &t) const { return Vec(x - t.x, y - t.y); }
	Vec operator * (const double &t) const { return Vec(x * t, y * t); }
	Vec turnLeft() const { return Vec(-y, x); }
	Vec rotate(double a) { return Vec(cos(a) * x - sin(a) * y, cos(a) * y + sin(a) * x); }
	double operator * (const Vec &t) const { return x * t.x + y * t.y; }
	double operator ^ (const Vec &t) const { return x * t.y - y * t.x; }
	double len2() { return x * x + y * y; }
	double len() { return sqrt(len2()); }
	double ang() const { return atan2(y, x); }
	bool operator < (const Vec &t) const { return ang() < t.ang(); }
} que[maxn], poly[maxn];
int q, n;

bool inPoly(Vec p) {
	bool ans = 0;
	rep(i, 0, n - 1) {
		Vec a = poly[i] - p, b = poly[nxt(i)] - p;
		if(dcmp(a.y) * dcmp(b.y) >= 0) continue;
		if(a.y > b.y) swap(a, b);
		if(dcmp(a ^ b) <= 0) continue;
		ans ^= 1;
	}
	return ans;
}

namespace SegCir {
	int tot;
	Vec crs[2];
	void intersect(double r, Vec p1, Vec p2) {
		Vec dir = (p1 - p2).turnLeft();
		double a = dir.x, b = dir.y, c = -a * p1.x - b * p1.y;
		double delta = 4.0 * b * b * (a * a * r * r - c * c + b * b * r * r);
		tot = 0;
		if(dcmp(delta) < 0) return ;
		if(delta < 0) delta = 0;
		double x, y;
		tot = 2;
		x = (-2.0 * a * c + sqrt(delta)) / (2.0 * (a * a + b * b));
		if(!dcmp(b)) y = sqrt(r * r - x * x); else y = (-a * x - c) / b;
		crs[0] = Vec(x, y);
		x = (-2.0 * a * c - sqrt(delta)) / (2.0 * (a * a + b * b));
		if(!dcmp(b)) y = -sqrt(r * r - x * x); else y = (-a * x - c) / b;
		crs[1] = Vec(x, y);
		if(!dcmp((crs[0] - crs[1]).len())) tot = 1;
		for(int i = 0; i < tot; i++)
			if(!dcmp((crs[i] - p2).len()) || dcmp((crs[i] - p1).len() + (crs[i] - p2).len() - (p1 - p2).len()))
				swap(crs[i], crs[tot-1]), i--, tot--;
		return ;
	}
}

Vec cs[maxn<<1];
int cnt;
double calc(double r) {
	if(!dcmp(r)) return inPoly(Vec(0, 0)) ? 1 : 0;
	cnt = 0;
	rep(i, 0, n - 1) {
		Vec a = poly[i], b = poly[nxt(i)];
		SegCir::intersect(r, a, b);
		rep(j, 0, SegCir::tot - 1) cs[cnt++] = SegCir::crs[j];
	}
	if(cnt <= 1) {
		double rang = pi / (rand() % 1000 + 1);
		Vec p(r * cos(rang), r * sin(rang));
		return inPoly(p) ? 1 : 0;
	}
	sort(cs, cs + cnt);
	double ans = 0;
	rep(i, 0, cnt - 1) {
		Vec a = cs[i], b = cs[(i+1)%cnt];
		double ang = a.ang(), bng = b.ang();
		if(bng < ang) bng += 2.0 * pi;
		double mng = (ang + bng) * .5;
		Vec mid(r * cos(mng), r * sin(mng));
		if(inPoly(mid)) ans += bng - ang;
	}
	return ans / (2.0 * pi);
}

int main() {
	srand((unsigned)time(NULL));
	
	int q = read(); n = read();
	double rang = pi / (rand() % 1000 + 1);
	rep(i, 1, q) {
		int x = read(), y = read();
		que[i] = Vec(x, y).rotate(rang);
	}
	rep(i, 0, n - 1) {
		int x = read(), y = read();
		poly[i] = Vec(x, y).rotate(rang);
	}
	
	double ans = 0;
	rep(i, 1, q) ans += calc(que[i].len());
	printf("%.5lf\n", ans);
	
	return 0;
}