動態幾何 史蒂芙的泡泡

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

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

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

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

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

將一個泡泡視為一個簡單多邊形 $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，便可以做到持久化的動態幾何操作，這便是一開始想要的目標。

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

參考解法

#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{
#defineKDMIN(a, b, c) {a[0] = min(b[0], c[0]), a[1] = min(b[1], c[1]);}
#defineKDMAX(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){
#definepushstk(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;
}