UOJ243【UR #16】破坏导蛋【计算几何,分块】
给定平面上 \(n\) 个点,\(m\) 次询问三角形内(包括边界上)有多少个点。
\(n\le 3\cdot 10^4\),\(m\le 10^5\)。
说到三角形数点,容斥一下可以变成角内数点,但这并没有什么用。
正解就是最暴力的半平面交:对于三个边界,分别算出对应范围的 bitset,然后 and 起来求 count。
跟 UOJ553 一样的套路,\(O(n^2)\) 维护按 \(y-kx\)(截距)排序的结果,所以还是分块。
设块大小是 \(B\),则复杂度为 \(O(nB\log n+\frac{nm\log n}B+\frac{nm}w)\),取 \(B=O(\sqrt m)\) 就得到复杂度是 \(O(n\sqrt m\log 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\n", ans[i]);
}