動態幾何 史蒂芙的泡泡 (解法 1)

contents

  1. 1. 題目描述
  2. 2. 解法
  3. 3. 參考解法

題目描述請至 Zerojudge e021: 史蒂芙的泡泡 查閱詳細內容

題目描述

在處理完數以百計的政事後,受盡折磨的史蒂芙,打算回家好好地休息。 拖著疲倦的身軀,再也無法再容納任何一點複雜計算。從王宮走回寢居的路上, 發現身邊所見的事物都不再圓滑,看起來就像是粗糙的幾何多邊形構成的一切。

打算享受著泡泡浴的史蒂芙,看著眼前的多邊形泡泡,失去原本應有的色澤,那透涼的心境更蒙上了一層灰影

「為什麼是我呢?」感嘆道

伸出手戳著眼前的泡泡,卻飄了過去

「區區的泡泡也跟我作對,嗚嗚」

將一個泡泡視為一個簡單多邊形 $A$,方便起見用一個序列 $a_0, a_1, ..., a_{n-1}$ 表示多邊形 $A$ 的每一個頂點,則會有 $n$ 個線段 $\overline{a_0 a_1}, \overline{a_1 a_2}, \cdots, \overline{a_{n-1} a_0}$

解法

從能找到的論文「A Unified Approach to Dynamic Point Location, Ray Shooting, and Shortest Paths in Planar Maps」 得知操作更新 $\mathcal{O}(\log^3 n)$,詢問操作 $\mathcal{O}(\log n)$,這一個實作難度估計沒辦法在 10K 限制下完成 (代碼上傳長度上限)。

從動態梯形剖分開始,複雜度就相當高了,而論文後要有類似輕重鏈剖分的概念,將對偶圖產生的節點以輕重邊保存,接著再維護雙線輪廓,讓相連的連續重邊可以保存左右兩側的輪廓,… 資訊量龐大,想到一些可能未提及的邊界案例,實作相當困難。而從其他題目的經驗中,當設計到 $\mathcal{O}(\log^2 n)$ 的操作時,由於操作常數大,運行時間可能無法匹敵 $\mathcal{O}(\sqrt{n})$ 的算法。

那麼有沒有其他的替代方案,只需要跑得比暴力法還要快就行?特別小心,暴力法找一個點是否在多邊形內,需要 $\mathcal{O}(n)$ 的時間,且快取效果非常好,很容易做到 SIMD 的平行技術。

首先,解決維護多邊形點集序列,可以使用任何的平衡樹完成,若採用 splay tree 在均攤 $\text{Amortized} \; \mathcal{O}(\log n)$。使用 treap 也可以完成這類操作,還可以做到額外的持久化功能。接著,修改點的操作,可以轉換成替換一個點的下一個點,因此將點放在 KD tree 上,而維護下一個點相當於把每一個節點擴充成 BRH。

接著,當解決點是否在多邊形內時,可走訪 BRH 找射線穿過的交點個數,此時複雜度至少為 $\mathcal{O}(\sqrt{n} + k)$,其中 $k$ 為交點個數。我們可以額外維護多邊形的順逆時針來優化搜尋空間,最後一個接觸的線段方向,額外提供奇偶數來判斷該點是否在內側,然後把射線縮短到交點到詢問點之間。不斷地縮短搜尋空間下,大部分的情況為 $\mathcal{O}(\sqrt{n})$,拋開了原本存在的 $k$

由於是封閉的簡單多邊形,所以當邊的疏密程度不一時,KD tree 轉換 BRH 會造成額外的負擔,bounding box 無法有效提供分割空間的功能。那麼在實際搜尋情況,額外走訪的節點與平面分布有關,這部分就不好分析。如果有更好的想法,歡迎分享。

現在實作動態 KD tree 必須有替罪羊樹 Scapegoat tree 的概念,再掛上 BRH 的想法,提供了相對有效率的動態方法來查詢點是否在多邊形內部。若把維護點序列的 splay tree 換成 treap,便可以做到持久化的動態幾何操作,這便是一開始想要的目標。

如果梯形剖分也可以持久化?暫時還不想去思考,那思考的容器要多大呢?

參考解法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
#include <bits/stdc++.h>
using namespace std;
struct Mesh {
static const int MAXN = 1e5 + 5;
int pt[MAXN][2];
void read(int n) {
for (int i = 0; i < n; i++)
scanf("%d %d", &pt[i][0], &pt[i][1]);
}
} mesh;
class SplayTree {
public:
struct Node {
Node *ch[2], *fa;
int size; int data;
Node() {
ch[0] = ch[1] = fa = NULL;
size = 1;
}
bool is_root() {
return fa->ch[0] != this && fa->ch[1] != this;
}
};
Node *root, *EMPTY;
void pushdown(Node *u) {}
void pushup(Node *u) {
if (u->ch[0] != EMPTY) pushdown(u->ch[0]);
if (u->ch[1] != EMPTY) pushdown(u->ch[1]);
u->size = 1 + u->ch[0]->size + u->ch[1]->size;
}
void setch(Node *p, Node *u, int i) {
if (p != EMPTY) p->ch[i] = u;
if (u != EMPTY) u->fa = p;
}
SplayTree() {
EMPTY = new Node();
EMPTY->fa = EMPTY->ch[0] = EMPTY->ch[1] = EMPTY;
EMPTY->size = 0;
}
void init() {
root = EMPTY;
}
Node* newNode() {
Node *u = new Node();
u->fa = u->ch[0] = u->ch[1] = EMPTY;
return u;
}
void rotate(Node *x) {
Node *y;
int d;
y = x->fa, d = y->ch[1] == x ? 1 : 0;
x->ch[d^1]->fa = y, y->ch[d] = x->ch[d^1];
x->ch[d^1] = y;
if (!y->is_root())
y->fa->ch[y->fa->ch[1] == y] = x;
x->fa = y->fa, y->fa = x;
pushup(y);
}
void deal(Node *x) {
if (!x->is_root()) deal(x->fa);
pushdown(x);
}
void splay(Node *x, Node *below) {
if (x == EMPTY) return ;
Node *y, *z;
deal(x);
while (!x->is_root() && x->fa != below) {
y = x->fa, z = y->fa;
if (!y->is_root() && y->fa != below) {
if (y->ch[0] == x ^ z->ch[0] == y)
rotate(x);
else
rotate(y);
}
rotate(x);
}
pushup(x);
if (x->fa == EMPTY) root = x;
}
Node* prevNode(Node *u) {
splay(u, EMPTY);
return maxNode(u->ch[0]);
}
Node* nextNode(Node *u) {
splay(u, EMPTY);
return minNode(u->ch[1]);
}
Node* minNode(Node *u) {
Node *p = u->fa;
for (; pushdown(u), u->ch[0] != EMPTY; u = u->ch[0]);
splay(u, p);
return u;
}
Node* maxNode(Node *u) {
Node *p = u->fa;
for (; pushdown(u), u->ch[1] != EMPTY; u = u->ch[1]);
splay(u, p);
return u;
}
Node* findPos(int pos) { // [0...]
for (Node *u = root; u != EMPTY;) {
pushdown(u);
int t = u->ch[0]->size;
if (t == pos) {
splay(u, EMPTY);
return u;
}
if (t > pos)
u = u->ch[0];
else
u = u->ch[1], pos -= t + 1;
}
return EMPTY;
}
tuple<int, int, int> insert(int data, int pos) { // make [pos] = data
Node *p, *q = findPos(pos);
Node *x = newNode(); x->data = data;
if (q == EMPTY) {
p = maxNode(root), splay(p, EMPTY);
setch(x, p, 0);
splay(x, EMPTY);
} else {
splay(q, EMPTY), p = q->ch[0];
setch(x, p, 0), setch(x, q, 1);
setch(q, EMPTY, 0);
splay(q, EMPTY);
p = prevNode(x);
}
if (p == EMPTY) p = maxNode(root);
if (q == EMPTY) q = minNode(root);
return make_tuple(p->data, data, q->data);
}
tuple<int, int, int> remove(int pos) {
Node *x = findPos(pos), *p, *q;
p = prevNode(x), q = nextNode(x);
if (p != EMPTY && q != EMPTY) {
setch(p, q, 1);
p->fa = EMPTY, splay(q, EMPTY);
} else if (p != EMPTY) {
p->fa = EMPTY, root = p;
} else {
q->fa = EMPTY, root = q;
}
int del = x->data;
free(x);
if (p == EMPTY) p = maxNode(root);
if (q == EMPTY) q = minNode(root);
return make_tuple(p->data, del, q->data);
}
int size() {
return root == EMPTY ? 0 : root->size;
}
} mlist;
static inline int log2int(int x) {
return 31 - __builtin_clz(x);
}
static inline int64_t h(int p, int q) {
return (int64_t) mesh.pt[p][0]*mesh.pt[q][1] - (int64_t) mesh.pt[p][1]*mesh.pt[q][0];
}
struct KDBRH {
static constexpr double ALPHA = 0.75;
static constexpr double LOG_ALPHA = log2(1.0 / ALPHA);
struct Pt {
int d[2];
Pt() {}
Pt(int xy[]):Pt(xy[0], xy[1]) {}
Pt(int x, int y) {d[0] = x, d[1] = y;}
bool operator==(const Pt &x) const {
return d[0] == x.d[0] && d[1] == x.d[1];
}
static Pt NaN() {return Pt(INT_MIN, INT_MIN);}
int isNaN() {return d[0] == INT_MIN;}
};
struct PtP {
Pt p, q;
PtP(Pt p, Pt q): p(p), q(q) {}
};
struct cmpAxis {
int k;
cmpAxis(int k): k(k) {}
bool operator() (const PtP &x, const PtP &y) const {
return x.p.d[k] < y.p.d[k];
}
};
struct BBox {
#define KDMIN(a, b, c) {a[0] = min(b[0], c[0]), a[1] = min(b[1], c[1]);}
#define KDMAX(a, b, c) {a[0] = max(b[0], c[0]), a[1] = max(b[1], c[1]);}
int l[2], r[2];
BBox() {}
BBox(int a[], int b[]) {
KDMIN(l, a, b); KDMAX(r, a, b);
}
void expand(Pt p) {
KDMIN(l, l, p.d); KDMAX(r, r, p.d);
}
void expand(BBox b) {
KDMIN(l, l, b.l); KDMAX(r, r, b.r);
}
inline int raycast(double x, double fx, double y) {
return l[1] <= y && y <= r[1] && r[0] >= x && l[0] <= fx;
}
static BBox init() {
BBox b; b.l[0] = b.l[1] = INT_MAX, b.r[0] = b.r[1] = INT_MIN;
return b;
}
};
struct Node {
Node *lson, *rson;
Pt pt, qt;
BBox box;
int size; int8_t used;
Node() {}
void init() {
lson = rson = NULL;
size = 1, used = 1;
pt = qt = Pt::NaN();
}
bool hasBox() { return box.l[0] <= box.r[0]; }
void pushup() {
size = used;
if (lson) size += lson->size;
if (rson) size += rson->size;
pushupBox();
}
void pushupBox() {
BBox t = BBox::init();
if (!qt.isNaN())
t.expand(pt), t.expand(qt);
if (lson && lson->hasBox())
t.expand(lson->box);
if (rson && rson->hasBox())
t.expand(rson->box);
box = t;
}
double interpolate(double y) {
if (pt.d[1] == qt.d[1]) return pt.d[0];
return pt.d[0] + (qt.d[0] - pt.d[0]) * (y - pt.d[1]) / (qt.d[1] - pt.d[1]);
}
};
Node *root;
vector<PtP> A;
int64_t area;
Node _mem[262144];
int gc[262144], gci, memi;
Node* newNode() {
Node *u = gci >= 0 ? &_mem[gc[gci--]] : &_mem[memi++];
u->init();
return u;
}
void freeNode(Node *u) {
gc[++gci] = u-_mem;
}
void init() {
root = NULL, area = 0;
gci = -1, memi = 0;
}
void insert(tuple<int, int, int> e) {
int p, q, r; tie(p, q, r) = e;
insert(root, 0, Pt(mesh.pt[q]), Pt(mesh.pt[p]), log2int(size()) / LOG_ALPHA);
changeNode(root, 0, Pt(mesh.pt[r]), Pt(mesh.pt[q]));
area += h(p, q) + h(q, r) - h(p, r);
}
void remove(tuple<int, int, int> e) {
int p, q, r; tie(p, q, r) = e;
remove(root, 0, Pt(mesh.pt[q]), log2int(size()) / LOG_ALPHA);
changeNode(root, 0, Pt(mesh.pt[r]), Pt(mesh.pt[p]));
area -= h(p, q) + h(q, r) - h(p, r);
}
int RAY_T;
double RAY_X;
vector<double> X;
int inside(double x, double y) {
if (area == 0) return 0;
X.clear(), RAY_X = 1e+12, RAY_T = -1;
raycast(root, x, y);
if (RAY_T < 0)
return X.size()&1;
int pass = (area > 0) == RAY_T;
for (auto &x : X)
pass += x <= RAY_X;
return pass&1;
}
int size() { return root == NULL ? 0 : root->size; }
inline int isbad(Node *u, Node *v) {
if (u == root) return 1;
int l = v ? v->size : 0;
l = max(l, u->size-u->used-l);
return l > u->size * ALPHA;
}
Node* build(int k, int l, int r) {
if (l > r) return NULL;
int mid = (l + r)>>1;
Node *ret = newNode();
sort(A.begin()+l, A.begin()+r+1, cmpAxis(k));
while (mid > l && A[mid].p.d[k] == A[mid-1].p.d[k])
mid--;
tie(ret->pt, ret->qt) = tie(A[mid].p, A[mid].q);
ret->lson = build(!k, l, mid-1);
ret->rson = build(!k, mid+1, r);
ret->pushup();
return ret;
}
void flatten(Node *u) {
if (!u) return ;
flatten(u->lson);
flatten(u->rson);
if (u->used) A.emplace_back(u->pt, u->qt);
freeNode(u);
}
void changeNode(Node *u, int k, Pt x, Pt qt) {
if (!u) return;
if (x == u->pt) {
u->qt = qt, u->pushupBox();
return;
}
changeNode(x.d[k] < u->pt.d[k] ? u->lson : u->rson, !k, x, qt);
u->pushupBox();
}
void rebuild(Node* &u, int k) {
A.clear(), A.reserve(u->size);
flatten(u);
u = build(k, 0, A.size()-1);
}
bool insert(Node* &u, int k, Pt x, Pt y, int d) {
if (!u) {
u = newNode(), u->pt = x, u->qt = y, u->pushup();
return d <= 0;
}
if (x == u->pt) {
u->used = 1, u->qt = y, u->pushup();
return d <= 0;
}
auto &v = x.d[k] < u->pt.d[k] ? u->lson : u->rson;
int t = insert(v, !k, x, y, d-1);
u->pushup();
if (t && !isbad(u, v))
return 1;
if (t) rebuild(u, k);
return 0;
}
bool remove(Node* &u, int k, Pt x, int d) {
if (!u)
return d <= 0;
if (x == u->pt) {
if (u->lson || u->rson)
u->used = 0, u->qt = Pt::NaN(), u->pushup();
else
freeNode(u), u = NULL;
return d <= 0;
}
auto &v = x.d[k] < u->pt.d[k] ? u->lson : u->rson;
int t = remove(v, !k, x, d-1);
u->pushup();
if (t && !isbad(u, v))
return 1;
if (t) rebuild(u, k);
return 0;
}
inline int cast(Node *u, double x, double y) {
if (u->qt.isNaN() || (u->pt.d[1] > y) == (u->qt.d[1] > y))
return 0;
double tx = u->interpolate(y);
if (tx <= x || tx > RAY_X)
return 0;
RAY_X = tx, RAY_T = u->pt.d[1] < u->qt.d[1];
X.emplace_back(tx);
return 1;
}
Node* stk[128];
void raycast(Node *u, double x, double y) {
#define pushstk(u) {*p++ = u;}
Node **p = stk;
pushstk(u);
while (p > stk) {
u = *--p;
if (!u || !u->size || !u->box.raycast(x, RAY_X, y))
continue;
cast(u, x, y);
pushstk(u->rson);
pushstk(u->lson);
}
}
} mbrh;
int main() {
int n, m, cmd, x, pos;
double px, py;
scanf("%d %d", &n, &m);
mesh.read(n);
mlist.init(), mbrh.init();
for (int i = 0; i < m; i++) {
scanf("%d", &cmd);
if (cmd == 1) {
scanf("%d %d", &x, &pos);
mbrh.insert(mlist.insert(x, pos));
} else if (cmd == 2) {
scanf("%d", &x);
mbrh.remove(mlist.remove(x));
} else {
scanf("%lf %lf", &px, &py);
puts(mbrh.inside(px, py) ? "1" : "0");
}
}
return 0;
}