c++ - 最近点对 O(nlogn) 算法 - C++ 实现中某些数据的问题

Question

我对找到最近点的分而治之算法有疑问。我在此页面上检查了 C++ 实现 https://www.geeksforgeeks.org/closest-pair-of-points-onlogn-implementation/但是这段代码有问题。它在大多数情况下工作正常，但对于某些数据，此实现返回其他结果，而不是暴力方法。例如，我们取十个点 (x, y):

(795 981)
(1905 4574)
(8891 665)
(6370 1396)
(93 8603)
(302 7099)
(326 5318)
(4493 3977)
(429 8687)
(9198 1558)

对于此数据，O(n log n) 算法返回 944.298 而不是暴力给出的 346.341。为什么会这样？

这正是 geeksforgeeks 对我的示例数据的实现:

#include <iostream>
#include <float.h>
#include <stdlib.h>
#include <math.h>
using namespace std;

struct Point
{
    int x, y;
};

int compareX(const void* a, const void* b)
{
    Point *p1 = (Point *)a,  *p2 = (Point *)b;
    return (p1->x - p2->x);
}

int compareY(const void* a, const void* b)
{
    Point *p1 = (Point *)a,   *p2 = (Point *)b;
    return (p1->y - p2->y);
}


float dist(Point p1, Point p2)
{
    return sqrt( (p1.x - p2.x)*(p1.x - p2.x) +
                 (p1.y - p2.y)*(p1.y - p2.y)
    );
}

float bruteForce(Point P[], int n)
{
    float min = FLT_MAX;
    for (int i = 0; i < n; ++i)
        for (int j = i+1; j < n; ++j)
            if (dist(P[i], P[j]) < min)
                min = dist(P[i], P[j]);
    return min;
}

float min(float x, float y)
{
    return (x < y)? x : y;
}


float stripClosest(Point strip[], int size, float d)
{
    float min = d;  // Initialize the minimum distance as d 

    for (int i = 0; i < size; ++i)
        for (int j = i+1; j < size && (strip[j].y - strip[i].y) < min; ++j)
            if (dist(strip[i],strip[j]) < min)
                min = dist(strip[i], strip[j]);

    return min;
}

float closestUtil(Point Px[], Point Py[], int n)
{
    // If there are 2 or 3 points, then use brute force 
    if (n <= 3)
        return bruteForce(Px, n);

    // Find the middle point 
    int mid = n/2;
    Point midPoint = Px[mid];


    Point Pyl[mid+1];   // y sorted points on left of vertical line 
    Point Pyr[n-mid-1];  // y sorted points on right of vertical line 
    int li = 0, ri = 0;  // indexes of left and right subarrays 
    for (int i = 0; i < n; i++)
    {
        if (Py[i].x <= midPoint.x)
            Pyl[li++] = Py[i];
        else
            Pyr[ri++] = Py[i];
    }

    float dl = closestUtil(Px, Pyl, mid);
    float dr = closestUtil(Px + mid, Pyr, n-mid);

    // Find the smaller of two distances 
    float d = min(dl, dr);

    Point strip[n];
    int j = 0;
    for (int i = 0; i < n; i++)
        if (abs(Py[i].x - midPoint.x) < d)
            strip[j] = Py[i], j++;

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

float closest(Point P[], int n)
{
    Point Px[n];
    Point Py[n];
    for (int i = 0; i < n; i++)
    {
        Px[i] = P[i];
        Py[i] = P[i];
    }

    qsort(Px, n, sizeof(Point), compareX);
    qsort(Py, n, sizeof(Point), compareY);

    // Use recursive function closestUtil() to find the smallest distance 
    return closestUtil(Px, Py, n);
}

// Driver program to test above functions 
int main()
{
    Point P[] = {{795, 981}, {1905, 4574}, {8891, 665}, {6370, 1396}, {93, 8603}, {302, 7099},
                 {326, 5318}, {4493, 3977}, {429, 8687}, {9198, 1558}};
    int n = sizeof(P) / sizeof(P[0]);
    cout << closest(P, n) << std::endl;
    cout << bruteForce(P, n) << std::endl;
    return 0;
} 

有人知道这里出了什么问题吗？我已经尝试修复它几个小时了，但我真的不明白为什么会出现这个问题。

Answer 1

由于 Pyl 和 Pyr 的大小分别为 mid+1 和 n-mid-1，下面两行

float dl = closestUtil(Px,     Pyl, mid  );
float dr = closestUtil(Px+mid, Pyr, n-mid);

应该改写如下:

float dl = closestUtil(Px,       Pyl, mid+1  );
float dr = closestUtil(Px+mid+1, Pyr, n-mid-1);

此外，正如this site的源代码中所注释的那样，在上面的代码中，假设所有 x 坐标都是不同的。例如，如果所有 x 坐标都相同，则 li 从 0 递增到 n 并且在 Pyl[li++] 处发生意外行为= Py[i]; 当 li >= mid+1 时。

---

顺便说一句，C++ 规范中根本不存在 VLA(Variable Vength Arrays)。由于数组 Px、Py、Pyl 和 Pyr 是在具有自动存储持续时间的堆栈上创建的，因此它们的大小应该是在编译时确定。但是包括 GNU 编译器在内的一些 C++ 编译器支持 VLA 作为编译器扩展，并允许声明具有动态长度的 C 风格数组。 (如何为 VLA 分配内存是特定于实现的。)但 C++ 通过 std::vector 提供动态数组功能，这可能使我们的代码更具可读性、可移植性和健壮性。