最近点対

概要

2 次元の点集合を入力とし，距離が最小の 2 点を出力する．点集合を $x$ 座標でソートしたあと，右半分と左半分で分割統治する．

計算量

$O(n \log n)$

使い方

両端 (半開区間) のポインタまたはイテレータを渡す．
事前に $x$ 座標昇順にソートしておくこと！
呼出し後は $y$ 座標昇順にソートされている．

実装

template <class iter>
std::tuple<ld, P, P> closest_pair(iter left, iter right) {
    int n = distance(left, right);

    if (n == 1) {
        return std::make_tuple(std::numeric_limits<ld>::max(), *left, *left);
    }

    if (n == 2) {
        if (left[0].imag() > left[1].imag()) swap(left[0], left[1]);
        return std::make_tuple(norm(left[0] - left[1]), left[0], left[1]);
    }

    iter middle = std::next(left, n / 2);
    ld x = middle->real();
    auto d = std::min(closest_pair(left, middle), closest_pair(middle, right));
    std::inplace_merge(left, middle, right,
                       [](const P &a, const P &b) { return a.imag() < b.imag(); });

    std::vector<iter> around;
    for (iter i = left; i != right; ++i) {
        ld dx = fabs(i->real() - x);
        ld &opt = std::get<0>(d);
        if (dx * dx >= opt) continue;
        for (auto j = around.rbegin(); j != around.rend(); ++j) {
            ld dx = i->real() - (**j).real();
            ld dy = i->imag() - (**j).imag();
            if (dy * dy >= opt) break;
            ld norm = dx * dx + dy * dy;
            if (opt > norm) {
                d = std::make_tuple(norm, *i, **j);
            }
        }
        around.push_back(i);
    }
    return d;
};

検証

AOJ CGL5A http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2061948#1
SPOJ CLOPPAIR #18124516
- long double -> WA, long long -> AC
- ちょうど6桁出力を忘れない

参考文献

蟻本