洛谷 P4023 - [CTSC2012] 极点统计（计算几何）

洛谷题面传送门

首先建出 \(S\) 的凸包，假设凸包上的点按顺时针排列的结果为 \(P_0,P_1,\dots,P_{k-1}\)。

不难发现往凸包加入一个点对凸包产生的贡献是将凸包变成其与一个三角形的并，具体来说，假设加入的点为 \(A\)，那么如果 \(A\) 在凸包内则凸包不变，否则假设其与 \(A\) 的两个切点为 \(P_x,P_y\)，那么凸包会变为其与 \(\triangle AP_xP_y\) 的并。

接下来考虑什么样的点是极点，显然如果扣掉凸包内的点，那么一个点不是极点当且仅当它不在任何一个除自己之外的 \(\triangle AP_xP_y\)。容易发现一个性质，就是如果我们把这些三角形它当作一个个环上的区间，那么如果一个区间 \([l_1,r_1]\) 完全被另一个区间 \([l_2,r_2]\) 严格包含（即 \(l_2<l_1<r_1<r_2\)），那么这个区间对应的点肯定不是极点，也就是说如果我们断环成链，那么所有极点对应的区间必然随着左端点的递增，右端点也是不降的，因此排个序拿个单调栈维护一下即可。

接下来就是实现的问题了。比较难实现的地方在于求一个点与凸包的切线，直接二分是不可取的，因为不存在单调性，正确的方法是，将凸包劈成上凸包和下凸包，首先我们可以对于凸包上一条线段 \(P_iP_{i+1}\)，求出其是否在三角形内，这可以用类似于半平面交中判断一个点在直线左边还有右边求出，这时候分类讨论，如果上凸包和下凸包上都有线段在三角形内（此时其要么都是最左边的线段，要么都是最右边的线段，只用进行 \(\mathcal O(1)\) 次 check），就直接在上下凸包上各二分一遍。否则其要么只在上凸包上，要么只在下凸包上，我们就现在上凸包上二分找到横坐标最大且 \(\le x_A\) 的点 \(P_i\)，显然如果 \(P_iP_{i+1}\) 不在三角形中就只能在下凸包上，否则我们往左往右各二分一遍找到对应的范围即可。

时间复杂度 \(n\log n\)。

由于有板子且实现比较垃圾，所以写得很长。不建议阅读。

const double EPS = 1e-8;
int sgn(double x) {return ((x < -EPS) ? -1 : ((x < EPS) ? 0 : 1));}
struct point {
	double x, y;
	point() {x = y = 0;}
	point(double _x, double _y) {x = _x; y = _y;}
	friend point operator + (const point &X, const point &Y) {return point(X.x + Y.x, X.y + Y.y);}
	friend point operator - (const point &X, const point &Y) {return point(X.x - Y.x, X.y - Y.y);}
	friend point operator * (const point &X, const double &Y) {return point(X.x * Y, X.y * Y);}
	friend double operator | (const point &X, const point &Y) {return X.x * Y.y - X.y * Y.x;}
	friend double operator ^ (const point &X, const point &Y) {return X.x * Y.x + X.y * Y.y;}
	friend bool operator == (const point &X, const point &Y) {return sgn(X.x - Y.x) == 0 && sgn(X.y - Y.y) == 0;}
	friend bool operator != (const point &X, const point &Y) {return sgn(X.x - Y.x) != 0 || sgn(X.y - Y.y) != 0;}
	friend int get_dir(point u, point v) {return sgn(u | v);} // v is on clockwise of u
	double operator ~ () const {return sqrt(x * x + y * y);}
	double operator ! () const {return atan2(y, x);}
};
bool cmp_ang(point u, point v) {
	auto up = [&](point u) {return sgn(u.y) > 0 || (sgn(u.y) == 0 && sgn(u.x) > 0);};
	bool U = up(u), V = up(v); return (U == V) ? (get_dir(u, v) > 0) : U;
}
bool cmp_xy(point u, point v) {return sgn(u.x - v.x) ? (sgn(u.x - v.x) < 0) : (sgn(u.y - v.y) < 0);}
struct line { // directed ray
	point A, B;
	line() {}
	line(point _A, point _B) {A = _A; B = _B;}
	int dir(point x) {return get_dir(B - A, x - A);} // 0: on line, 1: left of A -> B, -1: right of A -> B
	friend bool para(line x, line y) {return get_dir(x.B - x.A, y.B - y.A) == 0;}
	friend bool coin(line x, line y) {return para(x, y) && x.dir(y.A) == 0;}
	friend point inter(line x, line y) {
		if (para(x, y)) {fprintf(stderr, "bad line intersection!\n"); exit(-1);}
		double U = ((y.B - x.A) | (y.A - x.A)), V = ((y.B - x.B) | (y.A - x.B));
		return x.A + (x.B - x.A) * (U / (U - V));
	}
};
struct segment {
	point A, B;
	segment() {}
	segment(point _A, point _B) {A = _A; B = _B;}
	operator line() const {return line(A, B);}
	friend bool para(segment x, segment y) {return para(line(x), line(y));}
	bool on_seg(point x) {
		return line(*this).dir(x) == 0 && sgn(x.x - min(A.x, B.x)) >= 0 && sgn(x.x - max(A.x, B.x)) <= 0 &&
		sgn(x.y - min(A.y, B.y)) >= 0 && sgn(x.y - max(A.y, B.y)) <= 0;
	}
	friend bool has_inter(segment x, segment y) {
		if (para(x, y)) return 0;
		point p = inter(line(x), line(y));
		if (!x.on_seg(p) || !y.on_seg(p)) return 0;
		return 1;
	}
	friend point get_inter(segment x, segment y) {
		if (!has_inter(x, y)) {fprintf(stderr, "bad segment intersection!\n"); exit(-1);}
		return inter(line(x), line(y));
	}
};
double area(point x, point y, point z) {return fabs((z - x) | (y - x));}
struct convex {
	vector<point> v, up, dw;
	convex() {}
	void construct(vector<point> vec) {
		sort(vec.begin(), vec.end(), cmp_xy);
		v.clear(); up.clear(), dw.clear();
		vector<int> stk; stk.pb(0);
		for (int i = 1; i < vec.size(); i++) {
			while (stk.size() > 1 && get_dir(vec[i] - vec[stk.back()], vec[stk[stk.size() - 2]] - vec[stk.back()]) != 1) stk.ppb();
			stk.pb(i);
		}
		for (int i = 0; i < stk.size(); i++) v.pb(vec[stk[i]]), dw.pb(vec[stk[i]]);
		stk.clear(); stk.pb(0);
		for (int i = 1; i < vec.size(); i++) {
			while (stk.size() > 1 && get_dir(vec[i] - vec[stk.back()], vec[stk[stk.size() - 2]] - vec[stk.back()]) != -1) stk.ppb();
			stk.pb(i);
		}
		for (int i = (int)(stk.size()) - 2; i >= 1; i--) v.pb(vec[stk[i]]);
		for (int i = 0; i < stk.size(); i++) up.pb(vec[stk[i]]);
//		for (point p : v) printf("(%.10lf, %.10lf)\n", p.x, p.y);
	}
	convex(vector<point> vec) {construct(vec);}
	double poly_area() {
		double res = 0;
		for (int i = 1; i + 1 < v.size(); i++) res += area(v[0], v[i], v[i + 1]);
		return res;
	}
	bool in(point x) {
		if (get_dir(v[1] - v[0], x - v[0]) < 0) return 0;
		if (get_dir(v.back() - v[0], x - v[0]) > 0) return 0;
		if (v.size() == 2) return segment(v[0], v[1]).on_seg(x);
		int l = 2, r = v.size() - 2, p = 1;
		while (l <= r) {
			int mid = l + r >> 1;
			if (get_dir(v[mid] - v[0], x - v[0]) > 0) p = mid, l = mid + 1;
			else r = mid - 1;
		}
		return get_dir(v[p + 1] - v[p], x - v[p]) >= 0;
	}
} H;
const int MAXN = 1e5;
int n, m, vis[MAXN + 5]; point a[MAXN + 5];
double getang(point X, point Y, point Z) {return acos(((Y - X) ^ (Z - X)) / (~(Y - X)) / (~(Z - X)));}
bool in_tri(point X, point Y, point Z, point A) {return convex(vector<point>{X, Y, Z}).in(A);}
int main() {
//	freopen("extreme.in", "r", stdin);
//	freopen("extreme.out", "w", stdout);
	scanf("%d%d", &n, &m); vector<point> vec;
	for (int i = 1; i <= n; i++) {double x, y; scanf("%lf%lf", &x, &y); vec.pb(point(x, y));}
	for (int i = 1; i <= m; i++) scanf("%lf%lf", &a[i].x, &a[i].y);
	H.construct(vec);
	vector<tuple<int, int, int> > v;
	for (int i = 1; i <= m; i++) {
		if (H.in(a[i])) {vis[i] = 1; continue;}
		auto chk = [&](point x, point y) {return line(x, y).dir(a[i]) == -1;};
		int A = chk(H.dw[0], H.dw[1]), B = chk(H.dw[H.dw.size() - 2], H.dw.back());
		int C = chk(H.up.back(), H.up[H.up.size() - 2]), D = chk(H.up[1], H.up[0]);
		int pos1 = 0, pos2 = 0;
		if (A && D) {
			int l = 1, r = H.up.size() - 1, p = 0;
			while (l <= r) {
				int mid = l + r >> 1;
				if (chk(H.up[mid], H.up[mid - 1])) p = mid, l = mid + 1;
				else r = mid - 1;
			}
			pos1 = H.dw.size() - 1 + H.up.size() - 1 - p;
			l = 1, r = H.dw.size() - 1, p = 0;
			while (l <= r) {
				int mid = l + r >> 1;
				if (chk(H.dw[mid - 1], H.dw[mid])) p = mid, l = mid + 1;
				else r = mid - 1;
			}
			pos2 = p;
		} else if (B && C) {
			int l = 1, r = H.up.size() - 1, p = 0;
			while (l <= r) {
				int mid = l + r >> 1;
				if (chk(H.up[mid], H.up[mid - 1])) p = mid, r = mid - 1;
				else l = mid + 1;
			}
			pos2 = H.dw.size() - 1 + H.up.size() - 1 - (p - 1);
			l = 1, r = H.dw.size() - 1, p = 0;
			while (l <= r) {
				int mid = l + r >> 1;
				if (chk(H.dw[mid - 1], H.dw[mid])) p = mid, r = mid - 1;
				else l = mid + 1;
			}
			pos1 = p - 1;
		} else {
			int l = 0, r = H.dw.size() - 1, p = 0;
			while (l <= r) {
				int mid = l + r >> 1;
				if (H.dw[mid].x <= a[i].x) p = mid, l = mid + 1;
				else r = mid - 1;
			}
			if ((p == H.dw.size() - 1) ? B : chk(H.dw[p], H.dw[p + 1])) {
				l = 0, r = p - 1; int x = p, y = p;
				while (l <= r) {
					int mid = l + r >> 1;
					if (chk(H.dw[mid], H.dw[mid + 1])) x = mid, r = mid - 1;
					else l = mid + 1;
				}
				pos1 = x; l = p + 1; r = H.dw.size() - 1;
				while (l <= r) {
					int mid = l + r >> 1;
					if (chk(H.dw[mid - 1], H.dw[mid])) y = mid, l = mid + 1;
					else r = mid - 1;
				}
				pos2 = y;
			} else {
				l = 0, r = H.up.size() - 1; p = 0;
				while (l <= r) {
					int mid = l + r >> 1;
					if (H.up[mid].x <= a[i].x) p = mid, l = mid + 1;
					else r = mid - 1;
				}
				l = 0, r = p - 1; int x = p, y = p;
				while (l <= r) {
					int mid = l + r >> 1;
					if (chk(H.up[mid + 1], H.up[mid])) x = mid, r = mid - 1;
					else l = mid + 1;
				}
				pos2 = H.dw.size() - 1 + H.up.size() - 1 - x; l = p + 1; r = H.up.size() - 1;
				while (l <= r) {
					int mid = l + r >> 1;
					if (chk(H.up[mid], H.up[mid - 1])) y = mid, l = mid + 1;
					else r = mid - 1;
				}
				pos1 = H.dw.size() - 1 + H.up.size() - 1 - y;
			}
		}
//		printf("! %d %d %d\n", i, pos1, pos2);
		if (pos1 > pos2) v.pb(mt(pos1, pos2 + H.v.size(), i));
		else v.pb(mt(pos1, pos2, i)), v.pb(mt(pos1 + H.v.size(), pos2 + H.v.size(), i));
	}
	auto cmp = [&](tuple<int, int, int> x, tuple<int, int, int> y) {
		if (get<0>(x) != get<0>(y)) return get<0>(x) < get<0>(y);
		if (get<1>(x) != get<1>(y)) return get<1>(x) < get<1>(y);
		return getang(H.v[get<0>(y) % H.v.size()], H.v[(get<0>(x) + 1) % H.v.size()], a[get<2>(x)])
			 > getang(H.v[get<0>(y) % H.v.size()], H.v[(get<0>(y) + 1) % H.v.size()], a[get<2>(y)]);
	};
	sort(v.begin(), v.end(), cmp); vector<tuple<int, int, int> > stk;
	auto get_in = [&](int id, tuple<int, int, int> t) {
		return in_tri(H.v[get<0>(t) % H.v.size()], H.v[get<1>(t) % H.v.size()], a[get<2>(t)], a[id]);
	};
	for (auto t : v) {
		if (stk.empty() || get<2>(t) == get<2>(stk.back()) || !get_in(get<2>(t), stk.back())) {
			while (!stk.empty() && get_in(get<2>(stk.back()), t)) {
				if (get<2>(t) != get<2>(stk.back())) vis[get<2>(stk.back())] = 1;
				stk.ppb();
			}
			stk.pb(t);
		} else vis[get<2>(t)] = 1;
	}
	int res = 0; for (int i = 1; i <= m; i++) res += (!vis[i]);
	printf("%d\n", res);
	return 0;
}