您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> 關於C語言

分而治之演算法經典：最近點對問題 O(n log n) 完整解析

編輯：關於C語言

分而治之演算法經典：最近點對問題 O(n log n) 完整解析

最近點對問題給定平面上 n 點，找出歐氏距離最近的點對。暴力法 O(n²)，分而治之可達 **O(n log n)**。本文詳細剖析分割策略、合併步驟關鍵幾何性質 (每點最多檢查 7 個候選)、完整 C++ 實作與複雜度證明。

問題定義與暴力基準

輸入**：n 個點 (x_i, y_i)。
輸出**：距離最近的點對 (p, q) 及距離 d。
暴力法**：雙重迴圈計算所有 n(n-1)/2 對距離，O(n²)。

分而治之策略

分割**：依 x 座標中位數劃分垂直線 x = x_m，分成左半集合 A、右半集合 B，大小各約 n/2。
遞歸求解**：分別求 A 內最近點對 d_A、B 內最近點對 d_B。設 d = min(d_A, d_B)。
合併 (關鍵)**：檢查**跨越分割線**的點對是否有距離 < d。
- 只有 **x 座標距離分割線 < d** 的點可能構成更近點對 → 形成寬度 2d 的垂直帶狀區。
- 將帶區內點按 y 座標排序。
- 幾何引理**：對帶區內任意點 p，只需檢查其後 **最多 7 個** 點 (在 y 座標差 < d 範圍內)。
回傳三者中最小距離對應的點對。

幾何引理證明：為何最多檢查 7 點？

在 d × 2d 矩形中，任意兩點距離 ≥ d (否則違背 d 定義)。將矩形劃分為 8 個 d/2 × d/2 小正方形，每小正方形對角線 d/√2 < d，**最多容納 1 點**。點 p 所在小正方形以上下相鄰共 8 格，扣除自身格，**最多 7 格可能有點**，故最多檢查 7 點。

資料結構設計：雙陣列 + 歸併排序

X[]：按 x 排序的點陣列 (主陣列，遞歸分割用)。
Y[]：按 y 排序的點陣列 (輔助陣列，合併步驟用)。
Z[]：暫存陣列，用於歸併 Y 維持 y 序。
點結構體：{int id; double x, y;}，Y 中額外存 int idx 指向 X 中對應點。

完整 C++ 實作

struct Point {
    int id; double x, y;
    bool operator<=(const Point &o) const { return x <= o.x; } // for X sort
};
struct PointY {
    int idx; double x, y; // idx -> index in X
    bool operator<=(const PointY &o) const { return y <= o.y; } // for Y sort
};

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

void closest_pair_rec(Point X[], PointY Y[], PointY Z[], int l, int r, Point &a, Point &b, double &d) {
    if (r - l == 1) { // 2 points
        a = X[l]; b = X[r]; d = dist(X[l], X[r]); return;
    }
    if (r - l == 2) { // 3 points: brute force
        double d1 = dist(X[l], X[l+1]), d2 = dist(X[l+1], X[r]), d3 = dist(X[l], X[r]);
        if (d1 <= d2 && d1 <= d3) { a=X[l]; b=X[l+1]; d=d1; }
        else if (d2 <= d3) { a=X[l+1]; b=X[r]; d=d2; }
        else { a=X[l]; b=X[r]; d=d3; }
        return;
    }

    int m = (l + r) / 2;
    // Split Y into Z[l..m] (left) and Z[m+1..r] (right) by idx
    int f = l, g = m + 1;
    for (int i = l; i <= r; ++i)
        (Y[i].idx > m ? Z[g++] : Z[f++]) = Y[i];
    // Recurse: swap Y/Z roles
    Point al, bl, ar, br; double dl, dr;
    closest_pair_rec(X, Z, Y, l, m, al, bl, dl);
    closest_pair_rec(X, Z, Y, m+1, r, ar, br, dr);
    // Merge results
    if (dr < dl) { a = ar; b = br; d = dr; } else { a = al; b = bl; d = dl; }
    // Restore Y[l..r] sorted by y (merge Z back to Y)
    merge(Z + l, Z + m + 1, Z + m + 1, Z + r + 1, Y + l, [](auto &p, auto &q){ return p.y <= q.y; });

    // Build strip: points with |x - X[m].x| < d
    int k = l;
    for (int i = l; i <= r; ++i)
        if (fabs(Y[i].x - X[m].x) < d) Z[k++] = Y[i];
    // Check strip pairs (sorted by y)
    for (int i = l; i < k; ++i)
        for (int j = i+1; j < k && Z[j].y - Z[i].y < d; ++j) {
            double dp = dist(X[Z[i].idx], X[Z[j].idx]);
            if (dp < d) { d = dp; a = X[Z[i].idx]; b = X[Z[j].idx]; }
        }
}

bool closest_pair(Point X[], int n, Point &a, Point &b, double &d) {
    if (n < 2) return false;
    sort(X, X + n); // by x
    PointY *Y = new PointY[n], *Z = new PointY[n];
    for (int i = 0; i < n; ++i) { Y[i] = {i, X[i].x, X[i].y}; }
    sort(Y, Y + n); // by y
    closest_pair_rec(X, Y, Z, 0, n-1, a, b, d);
    delete[] Y; delete[] Z;
    return true;
}