UOJ243【UR #16】破壞導蛋【計算幾何,分塊】
阿新 • • 發佈:2021-06-24
給定平面上 \(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)\)
#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]); }