問題定義

給定二維平面上 n 個點，找出歐氏距離最近的一對點。暴力法 O(n²) 在大數據下不可接受，分治算法可達到 O(n log n)。

算法核心思想

1. 預排序：按 x 坐標排序，並保留按 y 排序的副本
2. 分治遞歸：
   li>以中位數 x 分割點集為左/右兩半 li>遞歸求左半、右半最近對距離 dl、dr li>d = min(dl, dr)
3. 合併關鍵：只需檢查 **中線帶寬 2d 內** 的點。幾何證明：帶內點按 y 排序後，每點只需與後續 **最多 7 個點** 比較（鴿巢原理）。
4. 合併時對帶內點按 y 排序，雙重循環檢查，更新全局最小距離。

關鍵數據結構

struct Point {
    double x, y;
    bool operator==(const Point& o) const { return x==o.x && y==o.y; }
};

using Points = std::vector<Point>;

// 輔助比較器
bool cmpX(const Point& a, const Point& b) { return a.x < b.x; }
bool cmpY(const Point& a, const Point& b) { return a.y < b.y; }

inline double dist(const Point& a, const Point& b) {
    double dx = a.x - b.x, dy = a.y - b.y;
    return std::sqrt(dx*dx + dy*dy);
}

線性時間選擇算法（中位數）

為保證 O(n log n) 最壞情況，需 **中位數的中位數** 算法 (BFPRT) 選取主元，避免快排退化。代碼包含：分組插入排序、遞歸選中位數、partition 分割。

核心遞歸函數

// 返回 [l, r] 區間最近對距離，結果存入 res[2]
double closestPair(Points& pts, int l, int r, Points& res) {
    int n = r - l + 1;
    if (n == 2) {
        res[0] = pts[l]; res[1] = pts[r];
        if (cmpX(pts[r], pts[l])) std::swap(pts[l], pts[r]);
        return dist(pts[l], pts[r]);
    }
    if (n == 3) {
        // 插入排序 3 元素
        for (int i = l+1; i <= r; ++i) {
            Point tmp = pts[i];
            int j = i-1;
            while (j >= l && cmpX(tmp, pts[j])) {
                pts[j+1] = pts[j]; --j;
            }
            pts[j+1] = tmp;
        }
        double d1 = dist(pts[l], pts[l+1]);
        double d2 = dist(pts[l+1], pts[l+2]);
        if (d1 < d2) { res[0]=pts[l]; res[1]=pts[l+1]; return d1; }
        else { res[0]=pts[l+1]; res[1]=pts[l+2]; return d2; }
    }

    int mid = (l + r) >> 1;
    Point median = select(pts, l, r, mid - l, cmpX); // 線性選擇中位數

    Points resL(2), resR(2);
    double dl = closestPair(pts, l, mid, resL);
    double dr = closestPair(pts, mid+1, r, resR);
    double d = std::min(dl, dr);
    Points& best = (d == dl) ? resL : resR;
    res[0] = best[0]; res[1] = best[1];

    // 收集中線帶內點
    Points strip;
    strip.reserve(r - l + 1);
    for (int i = mid+1; i <= r; ++i)
        if (pts[i].x - pts[mid].x < d) strip.push_back(pts[i]);
    for (int i = mid; i >= l; --i)
        if (pts[mid+1].x - pts[i].x < d) strip.push_back(pts[i]);

    // 按 y 排序
    std::sort(strip.begin(), strip.end(), cmpY);

    // 檢查每點後續最多 7 個點
    for (size_t i = 0; i < strip.size(); ++i) {
        for (size_t j = i+1; j < strip.size() && j <= i+7 && strip[j].y - strip[i].y < d; ++j) {
            double nd = dist(strip[i], strip[j]);
            if (nd < d) {
                d = nd;
                res[0] = strip[i]; res[1] = strip[j];
            }
        }
    }
    return d;
}

主程序與測試

int main() {
    int n; std::cin >> n;
    Points pts(n);
    for (int i = 0; i < n; ++i) std::cin >> pts[i].x >> pts[i].y;

    Points res(2);
    double ans = closestPair(pts, 0, n-1, res);
    std::cout << "Nearest: " << res[0] << " & " << res[1] << '\n';
    std::cout << "Distance: " << ans << '\n';
    return 0;
}

複雜度分析

• 預排序：O(n log n)
• 遞歸關係：T(n) = 2T(n/2) + O(n) → O(n log n)
• 空間：O(n) 輔助數組 + 遞歸棧 O(log n)

工程優化提示

1. 小規模 (≤32) 直接暴力，避免遞歸開銷
2. 遞歸傳遞按 y 排序的副本，合併時 std::merge 線性合併，省去重複排序
3. 距離比較用平方距離避免 sqrt，最終再開根
4. 坐標為整數時可用 64 位整數運算，完全避免浮點