最近接点対問題：分割統治法による O(n log n) 解法と実装

2次元平面上に n 個の点が与えられたとき、ユークリッド距離が最も小さい点の組み合わせ（最近接点対）を求める問題を扱います。愚直な全探索では O(n²) の計算量がかかりますが、分割統治法を用いることで O(n log n) で解くことが可能です。

問題の概要

点の集合 P = {p₁, p₂, ..., pₙ} が与えられたとき、以下を満たすペア (pᵢ, pⱼ) を見つけます。

min distance = √((xᵢ - xⱼ)² + (yᵢ - yⱼ)²)

入力と出力の例

入力:
(2, 3), (12, 30), (40, 50), (5, 1), (12, 10), (3, 4)

出力:
最小距離: 1.41421  (点 (2,3) と (3,4) の間の距離 √2 ≈ 1.41421)

アルゴリズムの流れ（分割統治法）

1. 事前ソート: 全点を x 座標順（xSorted）と y 座標順（ySorted）の2つの配列にソートしておく。
2. 分割: xSorted の中央値 midPoint を基準に、点集合を左右2グループに分割する。このとき ySorted も x 座標を基準に左右へ振り分け、ソート順を維持した ySortedLeft, ySortedRight を作る。
3. 再帰的解法: 左右それぞれで最近接距離 leftDist, rightDist を再帰的に求める。点の数が3以下なら全探索（O(1)）で直接計算する。
4. マージ（ストリップ領域のチェック）: 左右の最小距離 d = min(leftDist, rightDist) を求め、中央の垂直線から x 方向に距離 d 以内にある点だけを抽出し、strip 配列（y座標順）を作る。
5. ストリップ内探索: strip 内では、各点から y 座標が d 以内にある後続の点（最大7個まで）のみをチェックすれば十分であることを利用し、O(n) で最近接距離 stripDist を求める。
6. 結果: min(d, stripDist) を返す。

計算量の解析

• ソート: O(n log n)
• 再帰式: T(n) = 2T(n/2) + O(n) → O(n log n) (マスター定理)
• 全体: O(n log n)

疑似コード

1. 基底ケース（全探索）

function findMinDist(points[], n):
    minDist = ∞
    for i = 0 to n-1:
        for j = i+1 to n-1:
            if distance(points[i], points[j]) < minDist:
                minDist = distance(points[i], points[j])
    return minDist

2. ストリップ領域のチェック

function stripClose(strip[], size, d):
    minDist = d
    for i = 0 to size-1:
        for j = i+1 to size-1 and (strip[j].y - strip[i].y) < minDist:
            if distance(strip[i], strip[j]) < minDist:
                minDist = distance(strip[i], strip[j])
    return minDist

3. メインの再帰関数

function findClosest(xSorted[], ySorted[], n):
    if n <= 3:
        return findMinDist(xSorted, n)
    
    mid = n / 2
    midPoint = xSorted[mid]
    
    // ySorted を左右に振り分け（x座標で判定）
    ySortedLeft  = points in ySorted with x <= midPoint.x
    ySortedRight = points in ySorted with x > midPoint.x
    
    leftDist  = findClosest(xSorted[0..mid], ySortedLeft,  mid)
    rightDist = findClosest(xSorted[mid..n], ySortedRight, n-mid)
    
    d = min(leftDist, rightDist)
    
    // ストリップ作成
    strip = []
    for point in ySorted:
        if abs(point.x - midPoint.x) < d:
            strip.append(point)
    
    return min(d, stripClose(strip, strip.length, d))

C++ 実装例

標準ライブラリの sort と構造体を用いた実装です。

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <limits>

using namespace std;

struct Point {
    int x, y;
};

// x座標でソート用
bool cmpX(const Point& a, const Point& b) { return a.x < b.x; }
// y座標でソート用
bool cmpY(const Point& a, const Point& b) { return a.y < b.y; }

// ユークリッド距離
double dist(const Point& a, const Point& b) {
    return sqrt((double)(a.x - b.x)*(a.x - b.x) + (double)(a.y - b.y)*(a.y - b.y));
}

// 少数点数での全探索
double bruteForce(const vector<Point>& pts) {
    double minDist = numeric_limits<double>::max();
    for (size_t i = 0; i < pts.size(); ++i)
        for (size_t j = i + 1; j < pts.size(); ++j)
            minDist = min(minDist, dist(pts[i], pts[j]));
    return minDist;
}

// ストリップ内の最近接点対探索
double stripClosest(vector<Point>& strip, double d) {
    double minDist = d;
    // strip は y座標順でソート済み
    for (size_t i = 0; i < strip.size(); ++i) {
        for (size_t j = i + 1; j < strip.size() && (strip[j].y - strip[i].y) < minDist; ++j) {
            minDist = min(minDist, dist(strip[i], strip[j]));
        }
    }
    return minDist;
}

// 分割統治法の核心部分
double closestUtil(vector<Point>& xSorted, vector<Point>& ySorted) {
    int n = xSorted.size();
    if (n <= 3) return bruteForce(xSorted);

    int mid = n / 2;
    Point midPoint = xSorted[mid];

    // ySorted を左右に分割
    vector<Point> yLeft, yRight;
    for (const auto& p : ySorted) {
        if (p.x <= midPoint.x) yLeft.push_back(p);
        else yRight.push_back(p);
    }

    // xSorted も左右に分割（インデックスで管理してもよいが、ここではコピーで簡略化）
    vector<Point> xLeft(xSorted.begin(), xSorted.begin() + mid);
    vector<Point> xRight(xSorted.begin() + mid, xSorted.end());

    double dl = closestUtil(xLeft, yLeft);
    double dr = closestUtil(xRight, yRight);
    double d = min(dl, dr);

    // ストリップ作成
    vector<Point> strip;
    for (const auto& p : ySorted) {
        if (abs(p.x - midPoint.x) < d)
            strip.push_back(p);
    }

    return min(d, stripClosest(strip, d));
}

// エントリーポイント
double closestPair(vector<Point> points) {
    vector<Point> xSorted = points;
    vector<Point> ySorted = points;
    sort(xSorted.begin(), xSorted.end(), cmpX);
    sort(ySorted.begin(), ySorted.end(), cmpY);
    return closestUtil(xSorted, ySorted);
}

int main() {
    vector<Point> P = {{2, 3}, {12, 30}, {40, 50}, {5, 1}, {12, 10}, {3, 4}};
    cout << "The minimum distance is " << closestPair(P) << endl;
    return 0;
}

実行結果

The minimum distance is 1.41421

ポイントまとめ

• 事前に x, y それぞれでソートしておくことで、分割・マージステップを O(n) で行える。
• ストリップ領域では、幾何学的性質（半径 d の円に最大6個までしか点が入らない）により、内側のループが定数回（最大7回）で済むため全体が O(n) になる。
• 実装時は配列のコピーを避け、インデックスやイテレータで範囲を渡すとより高速・省メモリになる。