【APIO2018】选圆圈（平面分块 | CDQ分治 | KDT）

Description

给定平面上的 (n) 个圆，用三个参数 ((x, y, R)) 表示圆心坐标和半径。

每次选取最大的一个尚未被删除的圆删除，并同时删除所有与其相切或相交的圆。

最后输出每个圆分别是被那个圆所删除的。

Hint

• (1le nle 3 imes 10^5)
• (0le |x|, |y|, R le 10^9)

Solution 1

后来在 Codeforces 上找到的官方题解 Link here。如果对题解中某些说明无法理解可以参考上述内容。做法参考：Link here。

有一个非常简单的 (O(n^2)) 暴力，由于每次都要扫一遍所有圆所以复杂度爆炸。

我们尝试剪枝，缩小枚举的范围。

若当前最大圆的半径为 (R)，那么我们做一个分块操作：将整个平面 划分为一个个方格，方格的边长为 (2R)，刚好“框住”这个最大圆。

对于当前这个圆，我们只要搜索 其所在格子及其相邻的 即可（两个格子相邻定义为所在行的距离不超过 1，且列距离也不超过 1，换言之，一个格子所有相邻的就是周围一圈 8 个）。显然这样是不会漏记的，因为所有圆的半径都不超过方格大小，那么一定不会出现两个圆相交或相切，但却不在相邻两个方格内。

但我们总不能每次都搜怎么大的格子，最大圆的大小如果非常大而其他圆又很小的话这个剪枝没有丝毫用处。因此我们引入一个 重构机制：设现在考虑到的圆的半径为 (R^prime)，原来方格大小为 (L)。若 (R^prime < L)，那么 重构整个方格，并以 (2R^prime) 作为新的方格大小 (L)。

重构的复杂度会不会有问题？观察到，当半径不足方格大小的一半时才会重构，重构之后如果又来一次那又得一半。于是整个算法不会有超过 (O(log R)) 次重构，而一次重构的复杂度可以达到 (O(nlog n))（排序时重构，搜索时二分），所以总共最多也就两只 (log)，不会有问题。

可以证明每个圆被检查的次数为常数。因为对于一个圆，一个方格中不会有太多的大小相似的大圆对它进行检查，这些大圆会在之前就会被互相消除掉。

总复杂度为 (O(nlog nlog R))。

但其实重构是可以做到一次 (O(n)) 的，这样复杂度就是优秀的 (O(n(log n+log R)))。具体可以参考上述链接。

Code for Solution 1

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Problem : APIO2018 LOJ #2586 选圆圈
 */
#include <algorithm>
#include <iostream>
#include <vector>

using namespace std;
const int N = 3e5 + 5;
const int inf = 1e9;

int n;
struct Triad {
    int x, y, r;
    inline void read() {
        cin >> x >> y >> r;
        x += inf, y += inf;
    }
} circ[N];
int ans[N];

int cir_ord[N];
vector<Triad> f;
vector<int> g[N];
int bsiz;

inline long long sqr(int x) {
    return x * 1ll * x;
}
inline bool connect(Triad& a, Triad& b) {
    return sqr(a.x - b.x) + sqr(a.y - b.y) <= sqr(a.r + b.r);
}

typedef vector<Triad>::iterator vecIt;
inline void init_blocks(int L) {
    if (bsiz && bsiz / 2 < L) return;
    bsiz = L; f.clear();
    for (int i = 1; i <= n; i++) if (!ans[i])
        f.push_back(Triad{circ[i].x / bsiz, circ[i].y / bsiz, i});
    sort(f.begin(), f.end(), [](Triad& a, Triad& b) {
        return a.x != b.x ? a.x < b.x : a.y < b.y;
    });
    int cnt = 0;
    for (vecIt i = f.begin(), j; i != f.end(); i = j, ++cnt) {
        for (j = i; j != f.end(); j++)
            if (i->x != j->x || i->y != j->y) break;
        f[cnt] = *i, g[cnt].clear();
        for (vecIt k = i; k != j; k++) g[cnt].push_back(k->r);
    }
    f.resize(cnt);
}

void solve(int p) {
    if (ans[p]) return;
    init_blocks(circ[p].r * 2);

    int x = circ[p].x / bsiz;
    int y = circ[p].y / bsiz;
    for (int i = x - 1; i <= x + 1; i++)
        for (int j = y - 1; j <= y + 1; j++) {
            if (i < 0 || j < 0) continue;
            int c = lower_bound(f.begin(), f.end(), Triad{i, j, 0}, [](Triad& a, const Triad& b) {
                return a.x != b.x ? a.x < b.x : a.y < b.y;
            }) - f.begin();
            if (c == int(f.size()) || f[c].x != i || f[c].y != j)
                continue;

            vector<int> buf;
            for (auto k : g[c])
                if (connect(circ[k], circ[p])) ans[k] = p;
                else buf.push_back(k);
            g[c].swap(buf);
        }
}

signed main() {
    ios::sync_with_stdio(false);

    cin >> n;
    for (int i = 1; i <= n; i++)
        circ[i].read(), cir_ord[i] = i;
    sort(cir_ord + 1, cir_ord + 1 + n, [](int& a, int& b) {
        return circ[a].r != circ[b].r ? circ[a].r > circ[b].r : a < b;
    });

    for (int i = 1; i <= n; i++)
        solve(cir_ord[i]);
    for (int i = 1; i <= n; i++)
        cout << ans[i] << ' ';
    return cout << endl, 0;
}

Solution 2

虽然也是有理有据的 (log^2) 但跑的真的很慢。 做法参考：Link here。

对于一个圆 ((x, y, R))，考虑如何得到比它大的与之相交或相切，并且 不是因为其他圆而被带走的圆。

这些圆之间 不可能又有公共部分，否则它们之间会互相消除。

有一个比较显然的东西：这些圆一定在矩形 ((x-R, y-R) sim (x + R, y + R)) 的 边框 上。不会出现在内部是因为这些圆是比当前圆大的。

对于 (x) 方向，我们维护一下前面的扫描线，那么前面每个圆会有 (x+R, x-R) 两条。然后向后做，每次加入是将 (y) 方向的两条扫描线加进去，查询是查询 (y) 方向的 前驱后继 即可。这个用平衡树（set）显然可以搞。

然而发现这样只是求出一个圆的答案。现实相当于要多组询问，那么我们 CDQ 分治！

首先肯定是按 (R) 排序。然后先递归左侧，处理完左侧大的对右侧小的影响，然后递归右侧。如果是递归到重点 (l=r) 的话，发现这里还没有删除那就自己删掉，标记一下这个不是被带走的。

主要怎么做左侧对右侧的影响。首先对于左边的，我们取出所有 不是被带走的 圆的扫描线，作为“修改”；右侧则取出全部，作为“查询”。修改部分，形如 (x-R) 的是“插入”，(x+R) 是“删除”，分别表示加入可能的答案和移除这个已经没有的答案。查询部分显然两边都是要查询的。

最后得到了一堆扫描线，然后按坐标大小的顺序去做。插入时将这个圆的 (y) 坐标插入，删除时也将这个圆的 (y) 坐标删除。查询就是 (y) 坐标的前驱后继即可。注意答案对编号取 min。

这样的复杂度是 (O(nlog^2 n)) 的，然而被随机选择 KDT 吊起来打 QAQ

Code for Solution 2

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Problem : APIO2018 LOJ #2586 选圆圈
 */
#include <algorithm>
#include <cstring>
#include <iostream>
#include <set>
#include <vector>

using namespace std;
const int N = 3e5 + 5;

int n;
struct circle {
    int x, y, R, idx;
} C[N];
int ans[N];
int fin[N];

inline long long sqr(int x) {
    return x * 1ll * x;
}
inline bool connect(circle& a, circle& b) {
    return sqr(a.x - b.x) + sqr(a.y - b.y) <= sqr(a.R + b.R);
}

struct event {
    int time, cmd, pid;
};
bool self[N];

void cdqSolve(int l, int r) {
    if (l == r) {
        if (ans[l] > n) ans[l] = l, self[l] = 1;
        return;
    }
    int mid = (l + r) >> 1;
    cdqSolve(l, mid);
    
    set<pair<int, int> > s;
    set<pair<int, int> >::iterator it;
    vector<event> v;

    for (int i = l; i <= mid; i++) if (self[i]) {
        v.push_back(event{C[i].x - C[i].R, 1, i}); // ins
        v.push_back(event{C[i].x + C[i].R, 3, i}); // del
    }
    for (int i = mid + 1; i <= r; i++) {
        v.push_back(event{C[i].x - C[i].R, 2, i});
        v.push_back(event{C[i].x + C[i].R, 2, i}); // qry
    }
    sort(v.begin(), v.end(), [](event& a, event& b) {
        return a.time != b.time ? a.time < b.time : a.cmd < b.cmd;
    });
    for (auto i : v)
        if (i.cmd == 1) {
            s.insert({C[i.pid].y, i.pid});
        } else if (i.cmd == 3) {
            s.erase({C[i.pid].y, i.pid});
        } else {
            it = s.lower_bound({C[i.pid].y, 0});
            if (it != s.end()) {
                if (connect(C[i.pid], C[it->second]))
                    ans[i.pid] = min(ans[i.pid], it->second);
            }
            if (it != s.begin()) {
                if (connect(C[i.pid], C[(--it)->second]))
                    ans[i.pid] = min(ans[i.pid], it->second);
            }
        }

    v.clear();
    for (int i = l; i <= mid; i++) if (self[i]) {
        v.push_back(event{C[i].y - C[i].R, 1, i}); // ins
        v.push_back(event{C[i].y + C[i].R, 3, i}); // del
    }
    for (int i = mid + 1; i <= r; i++) {
        v.push_back(event{C[i].y - C[i].R, 2, i});
        v.push_back(event{C[i].y + C[i].R, 2, i}); // qry
    }
    sort(v.begin(), v.end(), [](event& a, event& b) {
        return a.time != b.time ? a.time < b.time : a.cmd < b.cmd;
    });
    for (auto i : v)
        if (i.cmd == 1) {
            s.insert({C[i.pid].x, i.pid});
        } else if (i.cmd == 3) {
            s.erase({C[i.pid].x, i.pid});
        } else {
            it = s.lower_bound({C[i.pid].x, 0});
            if (it != s.end()) {
                if (connect(C[i.pid], C[it->second]))
                    ans[i.pid] = min(ans[i.pid], it->second);
            }
            if (it != s.begin()) {
                if (connect(C[i.pid], C[(--it)->second]))
                    ans[i.pid] = min(ans[i.pid], it->second);
            }
        }
    
    cdqSolve(mid + 1, r);
}

signed main() {
    ios::sync_with_stdio(false);
    
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> C[i].x >> C[i].y >> C[i].R;
        C[i].idx = i;
    }

    sort(C + 1, C + 1 + n, [](circle& a, circle& b) {
        return a.R != b.R ? a.R > b.R : a.idx < b.idx;
    });
    memset(ans, 0x3f, sizeof(ans));
    cdqSolve(1, n);

    for (int i = 1; i <= n; i++)
        fin[C[i].idx] = C[ans[i]].idx;
    for (int i = 1; i <= n; i++)
        cout << fin[i] << ' ';
    return cout << endl, 0;
}

Solution 3

既然是“二维平面”，那么可以 KDT 暴力搞。

首先一个圆心为 ((x, y))，半径为 (R) 的圆，我们可以粗略地认为是一个矩形区域 ((x-R, y-R)sim (x+R, y+R))。

那么考虑用 KDT 维护这些矩形。对于一个圆，搜索可能与之有交集的区域即可。

但这样复杂度是 (O( ext{玄学})) 而且可能被卡。

于是发扬人类智慧，将所有的点 随机旋转 一个角度。于是就可以过 LOJ 数据了。

不过 Luogu 不需要旋转 qwq。

思维难度、代码难度和速度全方位碾压 Sol 1 & 2。

Code for Solution 3

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Problem : APIO2018 LOJ #2586 选圆圈
 */
#include <algorithm>
#include <cmath>
#include <iostream>

using namespace std;
const double eps = 1e-3;
const int N = 3e5 + 5;
const int K = 2;

int n;
struct circle {
    double p[K]; int r;
};
pair<circle, int> C[N];
int ans[N];

struct area {
    double max[K], min[K];
};

inline area trans(circle a) {
    return area{{a.p[0] + a.r, a.p[1] + a.r}, {a.p[0] - a.r, a.p[1] - a.r}};
}

struct node {
    int lc, rc;
    area range;
    circle val;
    int index, vaild;
} t[N];
int total = 0;

inline void pushup(int x) {
    t[x].range = trans(t[x].val);
    t[x].vaild = t[t[x].lc].vaild + t[t[x].rc].vaild + 1;
    for (int i = 0; i < K; i++) {
        if (t[x].lc) {
            t[x].range.max[i] = max(t[x].range.max[i], t[t[x].lc].range.max[i]);
            t[x].range.min[i] = min(t[x].range.min[i], t[t[x].lc].range.min[i]);
        }
        if (t[x].rc) {
            t[x].range.max[i] = max(t[x].range.max[i], t[t[x].rc].range.max[i]);
            t[x].range.min[i] = min(t[x].range.min[i], t[t[x].rc].range.min[i]);
        }
    }
}

int build(int l, int r, int d) {
    if (l > r) return 0;
    int mid = (l + r) >> 1, x = ++total;
    nth_element(C + l, C + mid, C + r + 1, [&](pair<circle, int> a, pair<circle, int> b) {
        return a.first.p[d] < b.first.p[d];
    });
    t[x].val = C[mid].first;
    t[x].index = C[mid].second;
    t[x].lc = build(l, mid - 1, d ^ 1);
    t[x].rc = build(mid + 1, r, d ^ 1);
    return pushup(x), x;
}

inline bool outside(circle& a, area& b) {
    double x = a.p[0], y = a.p[1]; int r = a.r;
    if (b.min[0] - (x + r) >= eps) return true;
    if (b.min[1] - (y + r) >= eps) return true;
    if ((x - r) - b.max[0] >= eps) return true;
    if ((y - r) - b.max[1] >= eps) return true;
    return false;
}
inline double sqr(double x) {
    return x * x;
}
inline bool connect(circle& a, circle& b) {
    return sqr(a.r + b.r) - (sqr(a.p[0] - b.p[0]) + sqr(a.p[1] - b.p[1])) >= -eps;
}
void select(int x, pair<circle, int>& c) {
    if (!x || outside(c.first, t[x].range)) return;
    if (!ans[t[x].index] && connect(c.first, t[x].val))
        ans[t[x].index] = c.second;
    if (t[t[x].lc].vaild) select(t[x].lc, c);
    if (t[t[x].rc].vaild) select(t[x].rc, c);
    t[x].vaild = t[t[x].lc].vaild + t[t[x].rc].vaild + (!ans[t[x].index] ? 1 : 0);
}

inline void rotate(double theta) {
    double SIN = sin(theta), COS = cos(theta);
    for (int i = 1; i <= n; i++) {
        double x = C[i].first.p[0];
        double y = C[i].first.p[1];
        C[i].first.p[1] = x * SIN + y * COS;
        C[i].first.p[0] = x * COS - y * SIN;
    }
}

signed main() {
    ios::sync_with_stdio(false);

    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> C[i].first.p[0] >> C[i].first.p[1] >> C[i].first.r;
        C[i].second = i;
    }
    rotate(2.3);

    int root = build(1, n, 0);

    sort(C + 1, C + 1 + n, [](pair<circle, int>& a, pair<circle, int>& b) {
        return a.first.r != b.first.r ? a.first.r > b.first.r : a.second < b.second;
    });
    for (int i = 1; i <= n; i++)
        if (!ans[C[i].second]) select(root, C[i]);
    
    for (int i = 1; i <= n; i++)
        cout << ans[i] << ' ';
    return cout << endl, 0;
}