UVa 1683 - In case of failure

contents

1. 1. Problem
2. 2. Sample Input
3. 3. Sample Output
4. 4. Solution
   1. 4.1. KD Tree
   2. 4.2. Delaunay

Problem

對平面上，每一個點找到最鄰近點的距離，輸出歐幾里得距離的平方即可。

Sample Input

2 
10 
17 41 
0 34 
24 19 
8 28 
14 12 
45 5 
27 31 
41 11 
42 45 
36 27 
15 
0 0 
1 2 
2 3 
3 2 
4 0 
8 4 
7 4 
6 3 
6 1 
8 0 
11 0 
12 2 
13 1 
14 2 
15 0

Sample Output

200 
100 
149 
100 
149 
52 
97 
52 
360 
97 
5 
2 
2 
2 
5 
1 
1 
2 
4 
5 
5 
2 
2 
2 
5

Solution

KD Tree

計算幾何中的資料結構 KD-tree，雖然不算完美，一開始先進行 random search，先找到基礎的剪枝條件，隨後掛上啟發式進行搜索。

KD-tree 的效能在 D 度空間，區間詢問 (回報區域中的所有點或求區域極值) 的效能約為 O(n^(1-1/d))。我不會估算找鄰近點對的效能。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <vector>
#include <queue>
using namespace std;
#define MAXN 262144
#define MAXD 2
#define INF (1LL<<60)
class kdTree {
public:
    struct PointD {
        long long d[MAXD];
        int label;
    };
    struct Node {
        Node *lson, *rson;
        int label;
    };
    struct cmp {
        bool operator()(const PointD* x, const PointD* y) const {
            if (x->d[sortDidx] != y->d[sortDidx])
                return x->d[sortDidx] < y->d[sortDidx];
            return x->label < y->label;
        }
    };
    struct pQcmp {
        bool operator() (pair<long long, int> &a, pair<long long, int> &b) const {
            if (a.first != b.first)	return a.first < b.first;
            return a.second < b.second;
        }
    };
    static int sortDidx;
    priority_queue<pair<long long, int>, vector< pair<long long, int> >, pQcmp> pQ; // <dist, label>
    Node buf[MAXN], *root;
    PointD pt[MAXN], *A[MAXN];
    int bufsize, K, qM, n;
    long long Q[MAXD], max_dist;
    void prepare(int n) {
        this->n = n;
        bufsize = 0;
        for (int i = 0; i < n; i++)
            A[i] = &pt[i];
        root = build(0, 0, n - 1);
    }
    Node* build(int k, int l, int r) {
        if(l > r)	return NULL;
        int m = (l + r)/2;
        Node *ret = &buf[bufsize++];
        sortDidx = k;
        nth_element(A + l, A + m, A + r + 1, cmp());
        ret->label = A[m]->label, ret->lson = ret->rson = NULL;
        if(l != r) {
            ret->lson = build((k+1)%K, l, m-1);
            ret->rson = build((k+1)%K, m+1, r);
        }
        return ret;
    }
    long long dist(int idx) {
        long long ret = 0;
        for(int i = 0; i < K; i++)
            ret += (long long)(pt[idx].d[i] - Q[i]) * (pt[idx].d[i] - Q[i]);
        return ret;
    }
    long long h_func(long long h[]) {
        long long ret = 0;
        for(int i = 0; i < K; i++)	ret += h[i];
        return ret;
    }
    void findNearest(Node *u, int k, long long h[]) {
        if(u == NULL || h_func(h) >= max_dist)
            return;
        long long d = dist(u->label);
        if (d <= max_dist) {
            pQ.push(make_pair(d, u->label));
            while (pQ.size() >= qM + 1)
                max_dist = pQ.top().first, pQ.pop();
        }
        long long old_hk = h[k];
        if (Q[k] <= pt[u->label].d[k]) {
            if (u->lson != NULL)
                findNearest(u->lson, (k+1)%K, h);
            if (u->rson != NULL) {
                h[k] = (pt[u->label].d[k] - Q[k]) * (pt[u->label].d[k] - Q[k]);
                findNearest(u->rson, (k+1)%K, h);
                h[k] = old_hk;
            }
        } else {
            if(u->rson != NULL)
                findNearest(u->rson, (k+1)%K, h);
            if (u->lson != NULL) {
                h[k] = (pt[u->label].d[k] - Q[k]) * (pt[u->label].d[k] - Q[k]);
                findNearest(u->lson, (k+1)%K, h);
                h[k] = old_hk;
            }
        }
    }
    void randomSearch(int ban) {
        max_dist = INF;
        long long d;
        for (int i = 0; i < 100; i++) {
            int x = rand()%n;
            if (x == ban)	continue;
            d = dist(x);
            if (d <= max_dist) {
                pQ.push(make_pair(d, x));
                while (pQ.size() >= qM + 1)
                    max_dist = pQ.top().first, pQ.pop();
            }
        }
    }
    void query(long long p[], int M, int ban) {
        for (int i = 0; i < K; i++)
            Q[i] = p[i];
        max_dist = INF, qM = M;
        long long h[MAXD] = {};
        randomSearch(ban);
        findNearest(root, 0, h);
        vector<int> ret;
        while (!pQ.empty())
            ret.push_back(pQ.top().second), pQ.pop();
        for (int i = (int) ret.size() - 1; i >= 0; i--) {
            if (ret[i] == ban)	continue;
            printf("%lld\n", dist(ret[i]));
            break;
        }
    }
};
int kdTree::sortDidx;

kdTree tree;
int main() {
    int n;
    long long p[MAXD];
    int testcase;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d", &n);
        tree.K = 2;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < 2; j++)
                scanf("%lld", &tree.pt[i].d[j]);
            tree.pt[i].label = i;
        }
        tree.prepare(n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < 2; j++)
                p[j] = tree.pt[i].d[j];
            tree.query(p, 2, i);
        }
    }
    return 0;
}
/*
 2
 10
 17 41
 0 34
 24 19
 8 28
 14 12
 45 5
 27 31
 41 11
 42 45
 36 27
 15
 0 0
 1 2
 2 3
 3 2
 4 0
 8 4
 7 4
 6 3
 6 1
 8 0
 11 0
 12 2 
 13 1 
 14 2 
 15 0
 */ 

Delaunay

關於 Delaunay triangulation 在此就不說明，由於是 Voronoi Diagram，就可以找到所有相鄰的點。

效能是 O(n log n)，而平面上的邊為線性 O(n)，但這種做法必須離線。

#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <set>
#include <map>
#include <queue>
#include <vector>
#include <string>
#include <iostream>
#include <assert.h>
#include <string.h>
#include <list>
using namespace std;

#define eps 1e-8
#define MAXN (1048576)
#define MAXV 262144
struct Point {
    double x, y;
    int id;
    Point(double a = 0, double b = 0, int c = -1):
    x(a), y(b), id(c) {}
    Point operator-(const Point &a) const {
        return Point(x - a.x, y - a.y);
    }
    Point operator+(const Point &a) const {
        return Point(x + a.x, y + a.y);
    }
    Point operator*(const double a) const {
        return Point(x * a, y * a);
    }
    Point operator/(const double a) const {
        return Point(x / a, y / a);
    }
    bool operator<(const Point &a) const {
        if (fabs(x - a.x) > eps)	return x < a.x;
        if (fabs(y - a.y) > eps)	return y < a.y;
        return false;
    }
    bool operator==(const Point &a) const {
        return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
    }
    bool operator!=(const Point &a) const {
        return !(fabs(x - a.x) < eps && fabs(y - a.y) < eps);
    }
    void read(int id = -1) {
        this->id = id;
    }
    double dist(Point b) {
        return hypot(x - b.x, y - b.y);
    }
    double dist2(Point b) {
        return (x - b.x) * (x - b.x) + (y - b.y) * (y - b.y);
    }
    void print() {
        printf("point (%lf, %lf)\n", x, y);
    }
};
struct Point3D {
    double x, y, z;
    Point3D(double a = 0, double b = 0, double c = 0):
    x(a), y(b), z(c) {}
    Point3D(Point p) {
        x = p.x, y = p.y, z = p.x * p.x + p.y * p.y;
    }
    Point3D operator-(const Point3D &a) const {
        return Point3D(x - a.x, y - a.y, z - a.z);
    }
    double dot(Point3D a) {
        return x * a.x + y * a.y + z * a.z;
    }
};
struct Edge {
    int id;
    list<Edge>::iterator twin;
    Edge(int id = 0) {
        this->id = id;
    }
};

int cmpZero(double v) {
    if (fabs(v) > eps)	return v > 0 ? 1 : -1;
    return 0;
}

double cross(Point o, Point a, Point b) {
    return (a.x-o.x)*(b.y-o.y) - (a.y-o.y)*(b.x-o.x);
}
Point3D cross(Point3D a, Point3D b) { // normal vector
    return Point3D(a.y * b.z - a.z * b.y
                   , -a.x * b.z + a.z * b.x
                   , a.x * b.y - a.y * b.x);
}
int inCircle(Point a, Point b, Point c, Point p) {
    if (cross(a, b, c) < 0)
        swap(b, c);
    Point3D a3(a), b3(b), c3(c), p3(p);
    //	printf("%lf %lf %lf\n", a3.x, a3.y, a3.z);
    //	printf("%lf %lf %lf\n", b3.x, b3.y, b3.z);
    //	printf("%lf %lf %lf\n", c3.x, c3.y, c3.z);
    //	printf("%lf %lf %lf\n", p3.x, p3.y, p3.z);
    b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
    Point3D f = cross(b3, c3); // normal vector
    return cmpZero(p3.dot(f)); // check same direction, in: < 0, on: = 0, out: > 0
}
int intersection(Point a, Point b, Point c, Point d) { // seg(a, b) and seg(c, d)
    return cmpZero(cross(a, c, b)) * cmpZero(cross(a, b, d)) > 0
    && cmpZero(cross(c, a, d)) * cmpZero(cross(c, d, b)) > 0;
}
class Delaunay {
public:
    list<Edge> head[MAXV]; // graph
    Point p[MAXV];
    int n, rename[MAXV];
    void init(int n, Point p[]) {
        for (int i = 0; i < n; i++)
            head[i].clear();
        for (int i = 0; i < n; i++)
            this->p[i] = p[i];
        sort(this->p, this->p + n);
        for (int i = 0; i < n; i++)
            rename[p[i].id] = i;
        this->n = n;
        divide(0, n - 1);
    }
    void addEdge(int u, int v) {
        head[u].push_front(Edge(v));
        head[v].push_front(Edge(u));
        head[u].begin()->twin = head[v].begin();
        head[v].begin()->twin = head[u].begin();
    }
    void divide(int l, int r) {
        if (r - l <= 1) { // #point <= 2
            for (int i = l; i <= r; i++)
                for (int j = i+1; j <= r; j++)
                    addEdge(i, j);
            return;
        }
        int mid = (l + r) /2;
        divide(l, mid);
        divide(mid + 1, r);
        
        list<Edge>::iterator it;
        int nowl = l, nowr = r;
        
        //		printf("divide %d %d\n", l, r);
        for (int update = 1; update; ) { // find left and right convex, lower common tangent
            update = 0;
            Point ptL = p[nowl], ptR = p[nowr];
            for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
                Point t = p[it->id];
                double v = cross(ptR, ptL, t);
                if (cmpZero(v) > 0 || (cmpZero(v) == 0 && ptR.dist2(t) < ptR.dist2(ptL))) {
                    nowl = it->id, update = 1;
                    break;
                }
            }
            if (update)	continue;
            for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
                Point t = p[it->id];
                double v = cross(ptL, ptR, t);
                if (cmpZero(v) < 0 || (cmpZero(v) == 0 && ptL.dist2(t) < ptL.dist2(ptR))) {
                    nowr = it->id, update = 1;
                    break;
                }
            }
        }
        
        addEdge(nowl, nowr); // add tangent
        //		printf("add base %d %d\n", nowl, nowr);
        for (int update = 1; true;) {
            update = 0;
            Point ptL = p[nowl], ptR = p[nowr];
            int ch = -1, side = 0;
            for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
                //				ptL.print(), ptR.print(), p[it->id].print();
                if (cmpZero(cross(ptL, ptR, p[it->id])) > 0
                    && (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0))
                    ch = it->id, side = -1;
                //				printf("test L %d %d %d\n", nowl, it->id, inCircle(ptL, ptR, p[ch], p[it->id]));
            }
            for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
                if (cmpZero(cross(ptR, p[it->id], ptL)) > 0
                    && (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0))
                    ch = it->id, side = 1;
                //				printf("test R %d %d %d\n", nowr, it->id, inCircle(ptL, ptR, p[ch], p[it->id]));
            }
            if (ch == -1)	break; // upper common tangent
            //			printf("choose %d %d\n", ch, side);
            if (side == -1) {
                for (it = head[nowl].begin(); it != head[nowl].end(); ) {
                    if (intersection(ptL, p[it->id], ptR, p[ch])) {
                        head[it->id].erase(it->twin);
                        head[nowl].erase(it++);
                    } else
                        it++;
                }
                nowl = ch;
                addEdge(nowl, nowr);
            } else {
                for (it = head[nowr].begin(); it != head[nowr].end(); ) {
                    if (intersection(ptR, p[it->id], ptL, p[ch])) {
                        head[it->id].erase(it->twin);
                        head[nowr].erase(it++);
                    } else
                        it++;
                }
                nowr = ch;
                addEdge(nowl, nowr);
            }
        }
    }
    
    vector< pair<int, int> > getEdge() {
        vector< pair<int, int> > ret;
        list<Edge>::iterator it;
        for (int i = 0; i < n; i++) {
            for (it = head[i].begin(); it != head[i].end(); it++) {
                if (it->id < i)
                    continue;
                //				printf("DG %d %d\n", i, it->id);
                ret.push_back(make_pair(p[i].id, p[it->id].id));
            }
        }
        return ret;
    }
} tool;
#define INF (1LL<<60)
Point p[MAXV];
long long ret[MAXV];
int main() {
    int testcase, n;
    long long x, y;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d", &n);
        for (int i = 0; i < n; i++) {
            scanf("%lld %lld", &x, &y);
            p[i] = Point(x, y, i);
        }
        tool.init(n, p);
        vector< pair<int, int> > DG = tool.getEdge();
        for (int i = 0; i < n; i++) {
            ret[i] = INF;
        }
        for (int i = 0; i < DG.size(); i++) {
            x = DG[i].first, y = DG[i].second;
            long long v = (long long) (p[x].x - p[y].x) * (long long) (p[x].x - p[y].x) +
            (long long) (p[x].y - p[y].y) * (long long) (p[x].y - p[y].y);
            ret[x] = min(ret[x], v);
            ret[y] = min(ret[y], v);
        }
        for (int i = 0; i < n; i++) {
            printf("%lld\n", ret[i]);
        }
    }
    return 0;
}
/*
 2
 10
 17 41
 0 34
 24 19
 8 28
 14 12
 45 5
 27 31
 41 11
 42 45
 36 27
 15
 0 0
 1 2
 2 3
 3 2
 4 0
 8 4
 7 4
 6 3
 6 1
 8 0
 11 0
 12 2
 13 1
 14 2
 15 0
 */