UOJ243【UR #16】破坏导蛋【计算几何，分块】

给定平面上 (n) 个点，(m) 次询问三角形内（包括边界上）有多少个点。

(nle 3cdot 10^4)，(mle 10^5)。

---

说到三角形数点，容斥一下可以变成角内数点，但这并没有什么用。

正解就是最暴力的半平面交：对于三个边界，分别算出对应范围的 bitset，然后 and 起来求 count。

跟 UOJ553 一样的套路，(O(n^2)) 维护按 (y-kx)（截距）排序的结果，所以还是分块。

设块大小是 (B)，则复杂度为 (O(nBlog n+frac{nmlog n}B+frac{nm}w))，取 (B=O(sqrt m)) 就得到复杂度是 (O(nsqrt mlog n+frac{nm}w))。虽然 bitset 部分复杂度看上去很大但瓶颈还是前者。

#include<bits/stdc++.h>
#define PB emplace_back
#define fi first
#define se second
#define MP make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> pii;
const int N = 30003, M = 100003, B = 200;
const double PI = acos(-1), eps = 1e-11;
template<typename T>
void rd(T &x){
	int ch = getchar(); x = 0; bool f = false;
	for(;ch < '0' || ch > '9';ch = getchar()) f |= ch == '-';
	for(;ch >= '0' && ch <= '9';ch = getchar()) x = x * 10 + ch - '0';
	if(f) x = -x;
}
struct Node {
	int x, y;
	Node(int _x = 0, int _y = 0): x(_x), y(_y){}
	double ang(){
		double res = atan2(y, x);
		if(res < 0) res += 2 * PI;
		return res;
	}
	Node operator - (const Node &o) const {return Node(x - o.x, y - o.y);}
	LL operator * (const Node &o) const {return (LL)x * o.y - (LL)y * o.x;}
} a[N], b[B], c[M][3];
void rd(Node &o){rd(o.x); rd(o.y);}
int n, m, ans[M], pos[B], rnk[B];
pair<double, int> q[M*3];
vector<pair<double, pii> > v;
bitset<B> res[M], bs[B];
void work(int n){
	sort(b, b+n, [&](const Node &p, const Node &q){return p.y == q.y ? p.x > q.x : p.y > q.y;});
	v.resize(0);
	for(int i = 0;i < n;++ i)
		for(int j = 0;j < n;++ j) if(i != j)
			v.PB(MP((b[i]-b[j]).ang(), MP(i, j)));
	sort(v.begin(), v.end(), [&](const pair<double, pii> &p, const pair<double, pii> &q){return fabs(p.fi - q.fi) <= eps ? p.se < q.se : p.fi < q.fi;});
	bs[0].reset(); bs[0].set(0);
	pos[0] = rnk[0] = 0;
	for(int i = 1;i < n;++ i){
		bs[i] = bs[i-1]; bs[i].set(i);
		pos[i] = rnk[i] = i;
	}
	int now = 0;
	for(int i = 0;i < m;++ i) res[i].set();
	for(int i = 0;i < m*3;++ i){
		while(now < v.size() && v[now].fi - q[i].fi <= eps){
			int x = v[now].se.fi, y = v[now].se.se, xx = rnk[x], yy = rnk[y];
			bs[xx].flip(x); bs[xx].flip(y);
			swap(rnk[x], rnk[y]); swap(pos[xx], pos[yy]); ++ now;
		} int id = q[i].se;
		Node p1 = c[id/3][id%3], p2 = c[id/3][(id+1)%3] - p1;
		int l = 0, r = n-1, md, rs = -1;
		while(l <= r){
			md = l+r>>1;
			if(p2*(b[pos[md]]-p1)>=0){rs = md; l = md+1;}
			else r = md-1;
		}
		if(~rs) res[id/3] &= bs[rs];
		else res[id/3].reset();
	}
	for(int i = 0;i < m;++ i)
		ans[i] += res[i].count();
}
int main(){
	rd(n); rd(m);
	for(int i = 0;i < n;++ i) rd(a[i]);
	for(int i = 0;i < m;++ i){
		for(int j = 0;j < 3;++ j) rd(c[i][j]);
		if((c[i][1]-c[i][0])*(c[i][2]-c[i][0]) < 0) swap(c[i][1], c[i][2]);
	}
	for(int i = 0;i < m;++ i){
		q[3*i] = MP((c[i][1]-c[i][0]).ang(), 3*i);
		q[3*i+1] = MP((c[i][2]-c[i][1]).ang(), 3*i+1);
		q[3*i+2] = MP((c[i][0]-c[i][2]).ang(), 3*i+2);
	}
	sort(q, q+3*m);
	for(int l = 0, r;l < n;l = r){
		r = min(l+B, n);
		for(int i = l;i < r;++ i) b[i-l] = a[i];
		work(r-l);
	}
	for(int i = 0;i < m;++ i) printf("%d
", ans[i]);
}