最近点对问题

编程之美2.11节。

也可参考：http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=113&page=show_problem&problem=1186

对于扩展问题2，可以先找到这些点的凸包(O(nlogn))，然后对凸包点求最大距离。对于凸包点，可以考虑按长轴平分四份，然后在最外侧的两份之间寻找最大点对。因为居内侧的两份中的最大距离无法超过长轴。当然，可以分割得更精细一点。这只是计算上复杂性的小差别了。

struct Point {
    double x, y;
};

double dist(const Point& p1, const Point& p2) {
    double dx = p1.x - p2.x, dy = p1.y - p2.y;
    return sqrt(dx*dx + dy*dy);
}

struct Compare {
    bool x;
    Compare(bool xx): x(xx) {}
    bool operator()(const Point& p, const Point& q) {
        return x? p.x < q.x: p.y < q.y;
    }
};

double closest(vector<Point>::iterator first, vector<Point>::iterator last) {
    typedef vector<Point>::iterator It;
    double r = MAX_DIST;
    if (last - first < 64) {
        for (It it = first; it != last; ++it) {
            for (It it2 = it + 1; it2 != last; ++it2) {
                double d = dist(*it, *it2);
                if (d < r) r = d;
            }
        }
        return r;
    }

    Point pmin, pmax;
    pmin = pmax = *first;
    for (It it = first + 1; it != last; ++it) {
        if (it->x < pmin.x) pmin.x = it->x;
        if (it->y < pmin.y) pmin.y = it->y;
        if (it->x > pmax.x) pmax.x = it->x;
        if (it->y > pmax.y) pmax.y = it->y;
    }

    bool px = pmax.y - pmin.y < pmax.x - pmin.x;
    It mid = first + (last - first)/2;
    nth_element(first, mid, last, Compare(px));

    double v = px? mid->x: mid->y;
    r = min(closest(first, mid), closest(mid, last));

    vector<Point> vec;
    for (It it = first; it != last; ++it)
        if (fabs(px? it->x - v: it->y - v) < r)
            vec.push_back(*it);
    sort(vec.begin(), vec.end(), Compare(!px));
    if (vec.end() - vec.begin() >= 8) {
        for (It it = vec.begin() + 8; it <= vec.end(); ++it) {
            double cr = closest(it - 8, it);
            if (cr < r) r = cr;
        }
    }
    else {
        double cr = closest(vec.begin(), vec.end());
        if (cr < r) r = cr;
    }

    return r;
}