土法炼钢兴趣小组的算法知识备份

最近点对与随机化几何算法

2025-07-15 · algorithms · #geometry #closest-pair #randomized #divide-and-conquer

目录

最近点对问题（Closest Pair Problem）是计算几何中最经典的问题之一：给定平面上 n 个点，找到欧氏距离最小的一对点。暴力枚举所有 C(n,2) 对的做法需要 O(n^2) 时间，而分治法可以将复杂度降至 O(n log n)——这一结果在 1975 年由 Shamos 和 Hoey 给出。更令人惊讶的是，Rabin 在 1976 年提出了一种随机化方法，可以在期望 O(n) 时间内解决这个问题。

本篇从暴力解法出发，逐步推导分治算法的正确性（特别是中间带最多检查 7 个邻居的证明），然后深入 Rabin 随机化算法的网格哈希思想，最后将视角拉高到随机增量构造（RIC）和 Clarkson-Shor 技术等通用框架，讨论高维最近邻的困难以及工程实践中的注意事项。

一、问题定义与暴力基线

1.1 形式化定义

给定平面上 n 个点的集合 P = {p_1, p_2, …, p_n}，其中 p_i = (x_i, y_i)。定义两点 p_i, p_j 之间的欧氏距离为：

d(p_i, p_j) = sqrt((x_i - x_j)^2 + (y_i - y_j)^2)

目标：找到一对点 (p_i, p_j)（i ≠ j）使得 d(p_i, p_j) 最小。

需要注意几个边界情形：

1.2 暴力算法

最直接的方法是枚举所有点对：

double brute_force(Point *pts, int n, int *idx1, int *idx2) {
    double min_dist = DBL_MAX;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            double d = dist_sq(pts[i], pts[j]);
            if (d < min_dist) {
                min_dist = d;
                *idx1 = i;
                *idx2 = j;
            }
        }
    }
    return min_dist;
}

时间复杂度 O(n^2)，空间复杂度 O(1)。当 n 在千万级别时，这个方法完全不可接受。

1.3 下界讨论

在代数决策树模型（algebraic decision tree model）下，最近点对问题的下界是 Omega(n log n)。这个下界可以通过规约元素唯一性问题（element uniqueness problem）来证明：如果我们能在 o(n log n) 时间内找到最近点对，就能在 o(n log n) 时间内判断一组数是否存在重复元素——而后者的下界是 Omega(n log n)。

然而，随机化算法可以绕过这个下界，因为代数决策树模型不允许随机性和哈希操作。Rabin 算法正是利用了这一点。

二、分治算法：O(n log n) 的经典方案

2.1 算法思路

分治法的核心思想如下：

  1. 预排序：将所有点按 x 坐标排序。
  2. 分割：用一条竖直线（通常取中位数点的 x 坐标）将点集分为左右两半，各包含 n/2 个点。
  3. 递归：分别递归求解左半部分和右半部分的最近点对，得到距离 delta_L 和 delta_R。
  4. 合并：令 delta = min(delta_L, delta_R)。最近点对可能跨越分割线，一个在左半部分、一个在右半部分。只需要在距离分割线不超过 delta 的”中间带”（strip）内寻找可能更近的点对。
  5. 中间带处理：将中间带内的点按 y 坐标排序，对每个点只需检查其后续的若干个邻居。

2.2 中间带最多检查 7 个邻居的证明

这是分治法的核心引理，也是整个算法正确性和复杂度分析的关键。

引理：在中间带内，对于按 y 坐标排序后的点序列，任意一个点 p 最多只需与其后续 7 个点比较。

证明：

考虑以分割线为中心、宽 2*delta 的竖直带。设 p 是中间带内的一个点，我们考虑 y 坐标在 [p.y, p.y + delta] 范围内且在中间带内的所有点。

这些点都落在一个 2*delta x delta 的矩形内。将这个矩形沿分割线分为两个 delta x delta 的正方形。

关键观察：在每个 delta x delta 的正方形内，任意两点之间的距离至少为 delta（因为同侧的最近距离已经是 delta_L 或 delta_R，而 delta = min(delta_L, delta_R)）。

在一个边长为 delta 的正方形内，能放下多少个两两距离至少为 delta 的点？答案是最多 4 个（将正方形等分为 4 个 delta/2 x delta/2 的小正方形，每个小正方形的对角线长度为 delta/sqrt(2) < delta，所以每个小正方形内最多一个点）。

因此：

这保证了中间带处理的时间复杂度是 O(n)，前提是中间带内的点已经按 y 排序。

2.3 递推关系与复杂度

算法的时间复杂度满足递推：

T(n) = 2T(n/2) + O(n)

其中 O(n) 来自中间带的处理（按 y 排序可以通过归并排序在递归过程中顺带完成，不需要额外的 O(n log n)）。

由主定理，T(n) = O(n log n)。

2.4 实现细节

分治法实现中有几个关键细节：

  1. 基础情形：当 n <= 3 时直接暴力计算。
  2. 归并排序与 y 排序的结合：在递归返回时，将左右两半按 y 坐标归并排序。这样中间带提取出来的点自然就是按 y 排序的。
  3. 浮点精度：比较距离时使用距离的平方，避免不必要的开方。
  4. 去重处理：如果有大量重合点，算法可以在找到距离为 0 时提前终止。

下面的示意图展示了分治过程中的中间带检查：

分治法中间带检查示意图

三、分治算法的完整 C 实现

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <string.h>

typedef struct {
    double x, y;
    int    id;       /* 原始编号，方便追踪 */
} Point;

typedef struct {
    int    i, j;     /* 最近点对的索引 */
    double dist_sq;  /* 距离的平方 */
} ClosestResult;

/* ---- 工具函数 ---- */

static inline double dist_sq(Point a, Point b) {
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return dx * dx + dy * dy;
}

static int cmp_x(const void *a, const void *b) {
    const Point *pa = (const Point *)a;
    const Point *pb = (const Point *)b;
    if (pa->x < pb->x) return -1;
    if (pa->x > pb->x) return  1;
    if (pa->y < pb->y) return -1;
    if (pa->y > pb->y) return  1;
    return 0;
}

static int cmp_y(const void *a, const void *b) {
    const Point *pa = (const Point *)a;
    const Point *pb = (const Point *)b;
    if (pa->y < pb->y) return -1;
    if (pa->y > pb->y) return  1;
    if (pa->x < pb->x) return -1;
    if (pa->x > pb->x) return  1;
    return 0;
}

/* ---- 暴力求解小规模 ---- */

static ClosestResult brute_force(Point *pts, int n) {
    ClosestResult res;
    res.dist_sq = DBL_MAX;
    res.i = res.j = -1;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            double d = dist_sq(pts[i], pts[j]);
            if (d < res.dist_sq) {
                res.dist_sq = d;
                res.i = pts[i].id;
                res.j = pts[j].id;
            }
        }
    }
    return res;
}

/* ---- 分治核心 ---- */

static ClosestResult closest_pair_rec(Point *px, int n, Point *tmp) {
    /* 基础情形：暴力求解 */
    if (n <= 3) {
        qsort(px, n, sizeof(Point), cmp_y);
        return brute_force(px, n);
    }

    int mid = n / 2;
    double mid_x = px[mid].x;

    /* 递归求解左右两半 */
    ClosestResult left  = closest_pair_rec(px, mid, tmp);
    ClosestResult right = closest_pair_rec(px + mid, n - mid, tmp);

    /* 取较小者 */
    ClosestResult best = (left.dist_sq <= right.dist_sq) ? left : right;

    /* 归并：按 y 坐标合并左右两半（此时两半各自已按 y 排序） */
    int i = 0, j = mid, k = 0;
    while (i < mid && j < n) {
        if (px[i].y <= px[j].y)
            tmp[k++] = px[i++];
        else
            tmp[k++] = px[j++];
    }
    while (i < mid) tmp[k++] = px[i++];
    while (j < n)   tmp[k++] = px[j++];
    memcpy(px, tmp, n * sizeof(Point));

    /* 提取中间带内的点 */
    int strip_cnt = 0;
    for (int s = 0; s < n; s++) {
        double dx = px[s].x - mid_x;
        if (dx * dx < best.dist_sq) {
            tmp[strip_cnt++] = px[s];
        }
    }

    /* 在中间带内检查——每个点最多看后续 7 个 */
    for (int s = 0; s < strip_cnt; s++) {
        for (int t = s + 1; t < strip_cnt && t <= s + 7; t++) {
            double dy = tmp[t].y - tmp[s].y;
            if (dy * dy >= best.dist_sq) break;
            double d = dist_sq(tmp[s], tmp[t]);
            if (d < best.dist_sq) {
                best.dist_sq = d;
                best.i = tmp[s].id;
                best.j = tmp[t].id;
            }
        }
    }

    return best;
}

/* ---- 对外接口 ---- */

ClosestResult closest_pair_dc(Point *pts, int n) {
    /* 复制一份，避免修改原数组 */
    Point *px  = (Point *)malloc(n * sizeof(Point));
    Point *tmp = (Point *)malloc(n * sizeof(Point));
    memcpy(px, pts, n * sizeof(Point));

    /* 按 x 排序 */
    qsort(px, n, sizeof(Point), cmp_x);

    ClosestResult res = closest_pair_rec(px, n, tmp);

    free(px);
    free(tmp);
    return res;
}

几个关键实现要点：

  1. 用 dist_sq 代替 dist 省去了开方运算，仅在最终输出时开方。
  2. 递归过程中顺带做归并排序，使得中间带的点自然按 y 排序，避免了额外的 O(n log n) 排序开销。
  3. tmp 数组在递归的不同层共享，但由于同一层不会同时使用同一段内存，所以是安全的——不过实现时需要小心。这里我们用了一个与 px 等长的临时数组，在归并和中间带提取中复用。
  4. 中间带检查用了 t <= s + 7 的循环上界和 dy * dy >= best.dist_sq 的提前退出条件，双重保障正确性。

四、Rabin 随机化算法：期望 O(n)

4.1 网格哈希的基本思想

Rabin 算法的核心观察是：如果我们已经知道了一个距离上界 delta，就可以将平面划分为边长为 delta 的网格。此时最近点对一定位于同一个网格单元或相邻的网格单元中。对于每个点，只需检查其所在网格和周围 8 个相邻网格中的点，就能找到比 delta 更近的点对（如果存在的话）。

将每个点映射到网格单元的操作是 O(1)（用哈希表存储），遍历所有点并检查邻居的总代价是 O(n)——因为每个点只检查常数个邻居格子。

问题在于：初始的 delta 从哪里来？

4.2 两轮采样策略

Rabin 算法使用两轮采样来获得一个好的初始 delta：

第一轮：从 n 个点中均匀随机采样 sqrt(n) 个点（不放回），暴力计算这些样本点的最近点对距离 delta_1。

构建网格：以 delta_1 为边长构建网格，将所有 n 个点插入网格。

第二轮：对于每个点，检查其所在网格和相邻 8 个网格中的所有点，更新最近距离。

算法伪代码：

RABIN-CLOSEST-PAIR(P):
    n = |P|
    随机打乱 P
    取前 m = ceil(sqrt(n)) 个点作为样本 S
    delta = BRUTE_FORCE(S)
    if delta == 0:
        return 0    /* 存在重合点 */
    构建边长为 delta 的网格 G
    将 P 中所有点插入 G
    for each point p in P:
        for each 相邻格子 c (包括自身所在格子):
            for each point q in c:
                if p != q:
                    delta = min(delta, dist(p, q))
    return delta

4.3 为什么期望线性：向后分析

证明 Rabin 算法期望 O(n) 的关键是向后分析（backward analysis）。

设真正的最近点对距离为 delta。我们需要证明：第一轮采样得到的 delta_1 不会比 delta 大太多，从而网格不会太粗糙，每个格子内的期望点数是 O(1)。

更精确地说，我们需要证明：当以 delta_1 为边长建网格时，格子中点的总数的期望是 O(n)。

向后分析的论证：

考虑将 n 个点按随机顺序逐一插入。设 delta_i 是前 i 个点的最近点对距离。当我们插入第 i 个点时，“当前的最近距离变小”这一事件等价于”第 i 个点参与了前 i 个点中的最近点对”。

对于 i 个随机排列的点中的最近点对，参与最近点对的两个点各有 1/i 的概率是最后插入的那个。因此最近距离因第 i 次插入而减小的概率最多为 2/i。

基于这个思想，可以推导出：

完整的分析需要更加细致的概率论证，但向后分析给出了直觉上的理解。

4.4 Rabin 算法的两阶段改进

上述基本版本有一个问题：第一轮采样可能得到很差的 delta_1（虽然概率很小），导致第二轮扫描实际运行时间很差。一种改进是使用两阶段策略：

  1. 第一阶段：随机选 sqrt(n) 个点，暴力求最近对距离 delta_1。
  2. 以 delta_1 建网格，插入所有点并扫描，得到更精确的 delta_2。
  3. 如果需要，以 delta_2 重建网格，再扫描一次。

实际上，一轮就足够了——期望复杂度已经是 O(n)。两阶段策略的意义在于减小常数因子和方差。

4.5 Rabin 算法的完整 C 实现

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <string.h>
#include <time.h>

typedef struct {
    double x, y;
    int    id;
} Point;

/* ---- 哈希网格 ---- */

typedef struct GridNode {
    int               pt_idx;   /* 指向 pts 数组的下标 */
    struct GridNode  *next;
} GridNode;

typedef struct {
    GridNode **buckets;
    int        capacity;
    double     cell_size;
} Grid;

/* 将 (gx, gy) 映射为哈希值 */
static unsigned int grid_hash(long long gx, long long gy, int cap) {
    /* 用一个简单的乘法哈希 */
    unsigned long long h = (unsigned long long)(gx * 1000000007LL + gy * 998244353LL);
    h ^= (h >> 16);
    h *= 0x45d9f3bULL;
    h ^= (h >> 16);
    return (unsigned int)(h % (unsigned long long)cap);
}

static Grid *grid_create(int capacity, double cell_size) {
    Grid *g = (Grid *)malloc(sizeof(Grid));
    g->capacity  = capacity;
    g->cell_size = cell_size;
    g->buckets   = (GridNode **)calloc(capacity, sizeof(GridNode *));
    return g;
}

static void grid_destroy(Grid *g) {
    for (int i = 0; i < g->capacity; i++) {
        GridNode *node = g->buckets[i];
        while (node) {
            GridNode *tmp = node->next;
            free(node);
            node = tmp;
        }
    }
    free(g->buckets);
    free(g);
}

static void grid_insert(Grid *g, Point *pts, int idx) {
    long long gx = (long long)floor(pts[idx].x / g->cell_size);
    long long gy = (long long)floor(pts[idx].y / g->cell_size);
    unsigned int h = grid_hash(gx, gy, g->capacity);
    GridNode *node = (GridNode *)malloc(sizeof(GridNode));
    node->pt_idx = idx;
    node->next   = g->buckets[h];
    g->buckets[h] = node;
}

/* 查询 (gx, gy) 格子中的所有点，与 pt 比较 */
static double grid_query_cell(Grid *g, Point *pts, int pt_idx,
                              long long gx, long long gy,
                              double best_sq, int *best_j)
{
    unsigned int h = grid_hash(gx, gy, g->capacity);
    for (GridNode *node = g->buckets[h]; node; node = node->next) {
        int j = node->pt_idx;
        if (j == pt_idx) continue;
        /* 需要确认这个点确实在 (gx, gy) 格子里（哈希冲突） */
        long long jgx = (long long)floor(pts[j].x / g->cell_size);
        long long jgy = (long long)floor(pts[j].y / g->cell_size);
        if (jgx != gx || jgy != gy) continue;

        double dx = pts[pt_idx].x - pts[j].x;
        double dy = pts[pt_idx].y - pts[j].y;
        double d  = dx * dx + dy * dy;
        if (d < best_sq) {
            best_sq = d;
            *best_j = j;
        }
    }
    return best_sq;
}

/* 查询点 pt_idx 周围 3x3 格子 */
static double grid_query_neighbors(Grid *g, Point *pts, int pt_idx,
                                   double best_sq, int *best_j)
{
    long long gx = (long long)floor(pts[pt_idx].x / g->cell_size);
    long long gy = (long long)floor(pts[pt_idx].y / g->cell_size);
    for (int dx = -1; dx <= 1; dx++) {
        for (int dy = -1; dy <= 1; dy++) {
            best_sq = grid_query_cell(g, pts, pt_idx,
                                      gx + dx, gy + dy,
                                      best_sq, best_j);
        }
    }
    return best_sq;
}

/* ---- Fisher-Yates 洗牌 ---- */

static void shuffle(int *arr, int n) {
    for (int i = n - 1; i > 0; i--) {
        int j = rand() % (i + 1);
        int tmp = arr[i];
        arr[i] = arr[j];
        arr[j] = tmp;
    }
}

/* ---- Rabin 随机化最近点对 ---- */

typedef struct {
    int    i, j;
    double dist_sq;
} RabinResult;

RabinResult closest_pair_rabin(Point *pts, int n) {
    RabinResult res;
    res.dist_sq = DBL_MAX;
    res.i = res.j = -1;

    if (n < 2) return res;

    srand((unsigned)time(NULL));

    /* 生成随机排列 */
    int *perm = (int *)malloc(n * sizeof(int));
    for (int i = 0; i < n; i++) perm[i] = i;
    shuffle(perm, n);

    /* 第一轮：取 sqrt(n) 个样本，暴力求最近对 */
    int m = (int)ceil(sqrt((double)n));
    if (m > n) m = n;

    double delta_sq = DBL_MAX;
    int bi = -1, bj = -1;
    for (int i = 0; i < m; i++) {
        for (int j = i + 1; j < m; j++) {
            double dx = pts[perm[i]].x - pts[perm[j]].x;
            double dy = pts[perm[i]].y - pts[perm[j]].y;
            double d  = dx * dx + dy * dy;
            if (d < delta_sq) {
                delta_sq = d;
                bi = perm[i];
                bj = perm[j];
            }
        }
    }

    if (delta_sq == 0.0) {
        res.dist_sq = 0.0;
        res.i = bi;
        res.j = bj;
        free(perm);
        return res;
    }

    double delta = sqrt(delta_sq);

    /* 第二轮：建网格，插入所有点，查询邻居 */
    int cap = 2 * n + 1;
    Grid *g = grid_create(cap, delta);
    for (int i = 0; i < n; i++) {
        grid_insert(g, pts, i);
    }

    res.dist_sq = delta_sq;
    res.i = bi;
    res.j = bj;

    for (int i = 0; i < n; i++) {
        int best_j = -1;
        double d = grid_query_neighbors(g, pts, i, res.dist_sq, &best_j);
        if (d < res.dist_sq) {
            res.dist_sq = d;
            res.i = i;
            res.j = best_j;
        }
    }

    grid_destroy(g);
    free(perm);
    return res;
}

五、随机增量构造（RIC）的一般范式

5.1 什么是 RIC

随机增量构造（Randomized Incremental Construction）是计算几何中一种通用的算法设计范式。其基本思想是：

  1. 将输入对象随机打乱。
  2. 逐一插入，每次插入后更新当前的”解”。
  3. 利用随机性保证每次更新的期望代价很小。

Rabin 的最近点对算法可以视为 RIC 的一个特例。更一般地，RIC 可以用于：

5.2 RIC 的一般模式

RIC(S):
    随机排列 S = {s_1, s_2, ..., s_n}
    初始化解 F_0 = 基于前 r 个对象的解（r 为常数）
    for i = r+1 to n:
        将 s_i 插入当前解 F_{i-1}
        if s_i 导致冲突（conflict）:
            更新 F_i（可能需要重建部分结构）
        else:
            F_i = F_{i-1}
    return F_n

5.3 冲突图（Conflict Graph）

RIC 的高效实现依赖于冲突图的维护。冲突图是一个二部图：

当新对象 s_i 被插入时：

  1. 通过冲突图找到所有与 s_i 冲突的组件（O(冲突数) 时间）。
  2. 删除这些组件，创建新组件。
  3. 更新冲突图（将未插入对象重新关联到新组件）。

如果每次更新的期望冲突数是常数，整个算法就是期望线性的。

5.4 最近点对的 RIC 视角

在最近点对问题中，RIC 可以这样理解：

关键问题：重建网格的代价是 O(i)（需要将前 i 个点重新插入新网格）。但期望重建次数是多少？

由向后分析，第 i 次插入导致最近距离变化的概率至多为 2/i。因此总期望重建代价为：

sum_{i=1}^{n} (2/i) * O(i) = sum_{i=1}^{n} O(1) = O(n)

这就是期望 O(n) 的来源。

六、Clarkson-Shor 技术：统一分析框架

6.1 动机

在许多计算几何问题中，我们需要分析”给定 n 个对象，某种几何结构的复杂度是多少”。Clarkson 和 Shor 在 1989 年提出了一种优雅的随机抽样框架，可以统一处理这类问题。

6.2 基本设定

设 H 是 n 个”半平面”（或超平面，更一般地是几何约束）的集合。对于 H 的一个子集 R，定义：

6.3 Clarkson-Shor 定理

定理：设 R 是从 H 中均匀随机选出的 r 个元素。设 f_0(R) 是 R 确定的 0 级顶点数量，f_{<=k}(H) 是 H 中级别不超过 k 的顶点数量。则：

E[f_0(R)] = Theta(f_{<=n/r}(H) * (r/n)^d)

这个定理的意义在于：通过随机抽样，我们可以从全集的组合复杂度推导出样本的期望复杂度，反之亦然。

6.4 应用举例

最近点对：

Delaunay 三角化：

线段求交：

6.5 与最近点对分析的联系

用 Clarkson-Shor 的语言重新表述 Rabin 算法的分析：

随机选取 r = sqrt(n) 个点作为样本 R。R 的最近距离 delta_R 满足：以 delta_R 为半径的圆盘内，全集 H 中期望有 O(n/r) = O(sqrt(n)) 个点。

因此，以 delta_R 为边长的网格中，每个格子期望只有 O(1) 个点。总共 n 个点，遍历邻居的总代价是 O(n)。

这种分析比直接的向后分析更加优雅，且能自然地推广到更复杂的问题。

七、高维最近邻的困难：维度诅咒

7.1 维度诅咒概述

上述 O(n log n) 的分治算法和 O(n) 的随机化算法都是针对二维平面的。当维度 d 增大时，情况会急剧恶化。

分治法的推广：

在 d 维空间中，分治法的中间带处理变为：每个点最多需要检查 O(C(d)) 个邻居，其中 C(d) 是一个与维度相关的常数，它随 d 指数增长。具体来说，中间带内的 delta x … x delta（d-1 维）的超矩形中能放下的两两距离至少为 delta 的点数，是关于 d 的指数函数。

分治法的递推关系变为：

T(n) = 2T(n/2) + O(n * C(d))

解为 O(C(d) * n log n)，其中 C(d) = 2^{O(d)}。

网格哈希的推广：

在 d 维空间中，每个格子有 3^d - 1 个邻居格子需要检查。当 d = 20 时，3^20 约等于 35 亿——这完全不可行。

7.2 近似最近邻

由于精确最近邻在高维下的困难，实际应用中常常转向近似最近邻（Approximate Nearest Neighbor, ANN）：

7.3 实际维度与内在维度

在实际数据中，虽然名义维度可能很高（如图像数据可能有数千维），但数据往往集中在一个低维流形上。这个”内在维度”（intrinsic dimension）通常远小于名义维度。许多数据结构的实际性能取决于内在维度而非名义维度。

倍增维度（doubling dimension）是一种常用的内在维度度量：一个度量空间的倍增维度是最小的 d，使得任何半径为 r 的球都可以被 2^d 个半径为 r/2 的球覆盖。在倍增维度为 d 的空间中，存在 2^{O(d)} * n log n 的最近邻算法。

八、k-d 树与球树

8.1 k-d 树

k-d 树是一种经典的空间划分数据结构，特别适合低维最近邻查询。

构造：

最近邻查询：

typedef struct KDNode {
    Point          pt;
    int            split_dim;
    struct KDNode *left;
    struct KDNode *right;
} KDNode;

/* 构建 k-d 树（二维情形） */
KDNode *kd_build(Point *pts, int n, int depth) {
    if (n <= 0) return NULL;

    int dim = depth % 2;  /* 0 = x, 1 = y */

    /* 选中位数作为根 */
    int mid = n / 2;
    /* 用 nth_element 的思路，这里简化为排序 */
    if (dim == 0)
        qsort(pts, n, sizeof(Point), cmp_x);
    else
        qsort(pts, n, sizeof(Point), cmp_y);

    KDNode *node = (KDNode *)malloc(sizeof(KDNode));
    node->pt        = pts[mid];
    node->split_dim = dim;
    node->left      = kd_build(pts, mid, depth + 1);
    node->right     = kd_build(pts + mid + 1, n - mid - 1, depth + 1);

    return node;
}

/* 最近邻查询 */
void kd_nearest(KDNode *root, Point query, double *best_sq, Point *best_pt) {
    if (!root) return;

    double d = dist_sq(root->pt, query);
    if (d < *best_sq) {
        *best_sq  = d;
        *best_pt  = root->pt;
    }

    int dim = root->split_dim;
    double diff = (dim == 0) ? query.x - root->pt.x : query.y - root->pt.y;

    KDNode *near_side  = (diff < 0) ? root->left  : root->right;
    KDNode *far_side   = (diff < 0) ? root->right : root->left;

    kd_nearest(near_side, query, best_sq, best_pt);

    /* 剪枝：分割面距离的平方 */
    if (diff * diff < *best_sq) {
        kd_nearest(far_side, query, best_sq, best_pt);
    }
}

8.2 球树（Ball Tree）

球树使用最小包围球来划分空间，而非轴对齐的超平面。

优点：

缺点：

8.3 k-d 树与球树的对比

方面 k-d 树 球树
分割方式 轴对齐超平面 最小包围球
构造时间 O(n log n) O(n log n)
查询（低维） O(log n) 期望 O(log n) 期望
查询（高维） 退化到 O(n) 比 k-d 树稍好
空间 O(n) O(n)
适用场景 低维、均匀分布 中维、非均匀分布
实现复杂度 简单 中等

8.4 与最近点对的关系

k-d 树和球树主要用于查询场景：给定一个查询点，找到数据集中离它最近的点。而最近点对问题是在整个数据集内部找距离最小的一对。

两者的联系：

九、工程实践与陷阱

9.1 工程陷阱汇总表

陷阱 症状 解决方案
距离开方 大量 sqrt 调用拖慢速度 全程使用距离平方比较，仅最终输出时开方
浮点精度 边界判断出错（如 dx*dx < delta_sq 误判） 加 epsilon 容差；或使用精确算术库
重合点 距离为 0 导致网格边长为 0、除零 特判：发现距离为 0 立即返回
哈希冲突 网格哈希中不同格子映射到同一桶 查询时验证格子坐标；使用更好的哈希函数
网格大小 哈希表容量不足导致链过长 容量设为 2n 以上，负载因子 < 0.5
整数溢出 格子坐标乘法溢出 使用 long long 或 double 计算格子坐标
随机种子 固定种子导致对特定输入退化 使用系统熵源（/dev/urandom）初始化
内存碎片 大量小块 malloc 导致缓存命中率低 预分配节点池（arena allocation）
递归深度 分治法在极端不平衡情况下栈溢出 限制递归深度；小规模切换暴力法
排序稳定性 不稳定排序影响相同 x/y 坐标点的处理 使用稳定排序或增加二级排序键
数据退化 所有点共线导致中间带包含所有点 分治法仍正确，但实际性能退化到 O(n log^2 n)
SIMD 优化 距离计算未向量化 使用 SOA 布局，利用 SSE/AVX 并行计算

9.2 性能优化策略

缓存友好的内存布局：

传统的 AOS（Array of Structures）布局：

/* AOS: 每个点的 x, y 相邻 */
Point pts[N];  /* {x0,y0}, {x1,y1}, ... */

SOA（Structure of Arrays）布局更利于 SIMD：

/* SOA: 所有 x 值连续，所有 y 值连续 */
double xs[N], ys[N];

当需要批量计算距离时，SOA 布局可以利用 SIMD 指令一次处理 4 个（SSE）或 8 个（AVX）点对的距离计算。

小规模切换阈值：

分治法在 n <= 某个阈值时切换到暴力法。最佳阈值取决于硬件，通常在 16-64 之间。过小则递归开销大，过大则暴力法的 O(n^2) 开始显现。

#define BRUTE_THRESHOLD 32

static ClosestResult closest_pair_rec(Point *px, int n, Point *tmp) {
    if (n <= BRUTE_THRESHOLD) {
        qsort(px, n, sizeof(Point), cmp_y);
        return brute_force(px, n);
    }
    /* ... 分治逻辑 ... */
}

并行化：

分治法天然适合并行化：左右两半的递归求解可以并行执行。使用 OpenMP 或线程池，可以在多核 CPU 上获得近线性加速。

/* OpenMP 并行版本（伪代码） */
#pragma omp parallel sections
{
    #pragma omp section
    left = closest_pair_rec(px, mid, tmp_l);
    #pragma omp section
    right = closest_pair_rec(px + mid, n - mid, tmp_r);
}

不过需要注意：并行化引入了额外的同步开销，只有在 n 足够大时才有意义。

十、算法对比与选择指南

10.1 三种算法的对比

维度 暴力法 分治法 Rabin 随机化
时间复杂度 O(n^2) O(n log n) O(n) 期望
最坏时间 O(n^2) O(n log n) O(n^2)
空间复杂度 O(1) O(n) O(n)
实现难度 简单 中等 较高（哈希表）
确定性 是 是 否
适用维度 任意 低维（d <= 10） 低维（d <= 5）
实际小规模性能 好（n < 100） 递归开销 哈希开销
常数因子 小 中等 较大
在线/离线 均可 离线 离线
并行友好 否 是 有限

10.2 选择建议

根据我的经验：

十一、完整测试驱动程序

下面给出一个完整的测试程序，同时运行分治法和 Rabin 算法并比较结果：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>

/* 假设上面的分治法和 Rabin 算法代码已经包含在内 */

/* ---- 测试生成器 ---- */

static Point *gen_random_points(int n, double range) {
    Point *pts = (Point *)malloc(n * sizeof(Point));
    for (int i = 0; i < n; i++) {
        pts[i].x  = ((double)rand() / RAND_MAX) * range;
        pts[i].y  = ((double)rand() / RAND_MAX) * range;
        pts[i].id = i;
    }
    return pts;
}

static Point *gen_collinear_points(int n) {
    Point *pts = (Point *)malloc(n * sizeof(Point));
    for (int i = 0; i < n; i++) {
        pts[i].x  = (double)i;
        pts[i].y  = 0.0;
        pts[i].id = i;
    }
    return pts;
}

static Point *gen_clustered_points(int n, int clusters) {
    Point *pts = (Point *)malloc(n * sizeof(Point));
    for (int i = 0; i < n; i++) {
        int c = i % clusters;
        double cx = c * 100.0;
        double cy = c * 100.0;
        pts[i].x  = cx + ((double)rand() / RAND_MAX) * 10.0;
        pts[i].y  = cy + ((double)rand() / RAND_MAX) * 10.0;
        pts[i].id = i;
    }
    return pts;
}

/* ---- 主测试 ---- */

int main(void) {
    srand((unsigned)time(NULL));

    int test_sizes[] = {10, 100, 1000, 10000};
    int num_tests = sizeof(test_sizes) / sizeof(test_sizes[0]);

    for (int t = 0; t < num_tests; t++) {
        int n = test_sizes[t];
        printf("=== 测试规模 n = %d ===\n", n);

        Point *pts = gen_random_points(n, 1000.0);

        /* 分治法 */
        clock_t start = clock();
        ClosestResult dc_res = closest_pair_dc(pts, n);
        clock_t end = clock();
        double dc_time = (double)(end - start) / CLOCKS_PER_SEC;

        /* Rabin 随机化 */
        start = clock();
        RabinResult rb_res = closest_pair_rabin(pts, n);
        end = clock();
        double rb_time = (double)(end - start) / CLOCKS_PER_SEC;

        printf("  分治法: dist = %.10f, 点 (%d, %d), 耗时 %.6f s\n",
               sqrt(dc_res.dist_sq), dc_res.i, dc_res.j, dc_time);
        printf("  Rabin:  dist = %.10f, 点 (%d, %d), 耗时 %.6f s\n",
               sqrt(rb_res.dist_sq), rb_res.i, rb_res.j, rb_time);

        /* 验证两种方法结果一致 */
        double diff = fabs(dc_res.dist_sq - rb_res.dist_sq);
        assert(diff < 1e-9 && "两种方法结果不一致！");
        printf("  验证通过。\n\n");

        free(pts);
    }

    /* 退化测试：共线点 */
    printf("=== 退化测试：共线点 ===\n");
    Point *collinear = gen_collinear_points(1000);
    ClosestResult col_res = closest_pair_dc(collinear, 1000);
    printf("  最近距离: %.10f (应为 1.0)\n", sqrt(col_res.dist_sq));
    assert(fabs(col_res.dist_sq - 1.0) < 1e-9);
    free(collinear);

    /* 退化测试：聚类点 */
    printf("=== 退化测试：聚类点 ===\n");
    Point *clustered = gen_clustered_points(5000, 10);
    ClosestResult cls_res = closest_pair_dc(clustered, 5000);
    printf("  最近距离: %.10f\n", sqrt(cls_res.dist_sq));
    free(clustered);

    printf("\n所有测试通过。\n");
    return 0;
}

十二、个人思考

12.1 算法选择的现实考量

在我多年的实践中，最近点对问题给我的最大感触是：理论最优不等于实践最优。

Rabin 算法在理论上是期望 O(n) 的，比分治法的 O(n log n) 更优。但在实际工程中，我很少推荐使用 Rabin 算法，原因如下：

  1. 常数因子大：哈希表的构建、维护和查询都有不小的常数开销。在 n < 10^6 时，分治法的实际运行时间往往更短。
  2. 最坏情况无保障：随机化算法有 O(n^2) 的最坏情况。虽然概率极小，但在安全关键（safety-critical）系统中不可接受。
  3. 实现复杂度高：哈希网格的正确实现需要处理大量边界情况（哈希冲突、整数溢出、内存管理等），远比分治法容易出错。
  4. 调试困难：随机化算法的非确定性使得 bug 难以复现。

分治法是我心中最近点对问题的”默认选择”。它确定、高效、实现清晰、易于调试、容易并行化。

12.2 随机化的真正价值

但这并不意味着随机化没有价值。恰恰相反，我认为随机化思想在计算几何中的价值是方法论层面的：

这些思想的价值远超最近点对这一个问题。理解了它们，你就掌握了一套分析和设计计算几何算法的通用工具。

12.3 高维的现实

维度诅咒是计算几何中最深刻的困难之一。当维度超过 10-15 时，几乎所有基于空间划分的精确算法都会退化到接近暴力枚举的性能。

这不是算法设计者的失败——这是问题本身的固有困难。高维空间的几何直觉完全不同于低维：所有点都近乎等距，“近邻”的概念变得模糊。

实际应用中的应对策略：

12.4 竞赛与面试

如果你是为了竞赛或面试准备，我的建议是：

  1. 必须掌握：分治法的思路和复杂度分析。中间带检查 7 个邻居的证明。能够手写分治法实现。
  2. 了解即可：Rabin 算法的思路、RIC 范式、Clarkson-Shor 技术。
  3. 不需要会：Rabin 算法的完整实现细节（竞赛中基本不会用到，面试中最多口述思路）。

最近点对是一个极好的”分治法教学样例”——它体现了分治法的三个关键步骤（分割、递归、合并），而合并步骤的巧妙处理（中间带的 O(n) 合并）是整个算法的精华。理解了这一点，你对分治法的理解就上了一个台阶。

参考文献

  1. Shamos, M. I., & Hoey, D. (1975). Closest-point problems. Proceedings of the 16th Annual Symposium on Foundations of Computer Science (FOCS), 151-162.
  2. Rabin, M. O. (1976). Probabilistic algorithms. Algorithms and Complexity: New Directions and Recent Results, 21-39.
  3. Clarkson, K. L., & Shor, P. W. (1989). Applications of random sampling in computational geometry, II. Discrete & Computational Geometry, 4(1), 387-421.
  4. Cormen, T. H., Leiserson, C. E., Rivest, R. L., & Stein, C. (2009). Introduction to Algorithms (3rd ed.). MIT Press. Chapter 33: Computational Geometry.
  5. de Berg, M., Cheong, O., van Kreveld, M., & Overmars, M. (2008). Computational Geometry: Algorithms and Applications (3rd ed.). Springer-Verlag.
  6. Mulmuley, K. (1994). Computational Geometry: An Introduction Through Randomized Algorithms. Prentice Hall.
  7. Har-Peled, S. (2011). Geometric Approximation Algorithms. American Mathematical Society.
  8. Indyk, P., & Motwani, R. (1998). Approximate nearest neighbors: towards removing the curse of dimensionality. Proceedings of the 30th Annual ACM Symposium on Theory of Computing (STOC), 604-613.
  9. Friedman, J. H., Bentley, J. L., & Finkel, R. A. (1977). An algorithm for finding best matches in logarithmic expected time. ACM Transactions on Mathematical Software, 3(3), 209-226.
  10. Kleinberg, J. M. (1997). Two algorithms for nearest-neighbor search in high dimensions. Proceedings of the 29th Annual ACM Symposium on Theory of Computing (STOC), 599-608.

算法系列导航：上一篇：点定位与梯形分解 | 下一篇：树形 DP

相关阅读：凸包算法 | 随机化算法

By .