|
2. It's shorter and quite a bit faster;
3. It's cleaner IMO.
<code snipped, needs more thinking>
<lang C>#include <stdio.h>
#include <stdlib.h>
#include <values.h>
#include <math.h>
typedef struct { double x, y; } point_t, *point;
/* note: even though l_list and r_list are used by each recursion of closest(),
they are always used in sequentially, no need to allocate them repeatedly */
point *l_list, *r_list;
inline double dist(point a, point b)
{
double dx = a->x - b->x, dy = a->y - b->y;
return dx * dx + dy * dy;
}
inline int cmp_dbl(double a, double b)
{
return a < b ? -1 : a > b ? 1 : 0;
}
double brute_force(point pts, int max_n, point *a, point *b)
{
int i, j;
double d, min_d = MAXDOUBLE;
for (i = 0; i < max_n; i++) {
for (j = i + 1; j < max_n; j++) {
d = dist(pts + i, pts + j);
if (d >= min_d ) continue;
*a = pts + i;
*b = pts + j;
min_d = d;
}
}
return min_d;
}
void sort_y(point *list, int n)
{
int cmp_y(const void *a, const void *b) {
return cmp_dbl( ((point)a)->y, ((point)b)->y );
}
qsort(list, n, sizeof(point), cmp_y);
}
double closest(point pts, int n, point *a, point *b)
{
int left, right, lsize, rsize, i;
point a1, b1;
double min_d, d, dsqrt, median;
/* problem small enough, don't divide, just conquer */
if (n <= 8) return brute_force(pts, n, a, b);
/* get left and right results */
left = right = n / 2;
min_d = closest(pts, left, a, b);
d = closest(pts + left, n - left, &a1, &b1);
if (min_d > d) {
min_d = d;
*a = a1;
*b = b1;
}
median = pts[left].x;
dsqrt = sqrt(min_d);
/* find points within +- min distance on X from the center line,
list their pointers */
for ( lsize = 0, left--;
median - pts[left].x < dsqrt && left;
l_list[lsize++] = pts + left--);
for ( rsize = 0;
pts[right].x - median < dsqrt && right < n;
r_list[rsize++] = pts + right++);
/* sort the pointers by y: don't touch the point data which is sorted by x */
sort_y(l_list, lsize);
sort_y(r_list, rsize);
/* climb up left and right list and compare distance */
for (left = right = 0; left < lsize; left ++) {
/* next point in left list */
median = l_list[left]->y;
/* climb up right list until the y is not too low */
while (right < rsize && median > dsqrt + r_list[right]->y)
right++;
if (right >= rsize) break;
for (i = right; i < rsize; i++) {
/* right y is too high, break */
if (r_list[i]->y > median + dsqrt) break;
if (min_d > (d = dist(r_list[i], l_list[left]))) {
*a = l_list[left];
*b = r_list[i];
min_d = d;
dsqrt = sqrt(min_d);
}
}
}
return min_d;
}
#define NP 1000000
int main(int argc, char **argv)
{
int i;
point a, b;
point pts = malloc(sizeof(point_t) * NP);
l_list = malloc(sizeof(point) * (NP / 2));
r_list = malloc(sizeof(point) * ((NP + 1) / 2));
for(i = 0; i < NP; i++) {
pts[i].x = 20 * (double) rand()/RAND_MAX;
pts[i].y = 20 * (double) rand()/RAND_MAX;
}
/*
printf("brute force: %g, ", sqrt(brute_force(pts, NP, &a, &b)));
printf("between (%f,%f) and (%f,%f)\n", a->x, a->y, a->x, a->y);
*/
int cmp_x(const void *a, const void *b) {
return cmp_dbl( ((point)a)->x, ((point)b)->x );
}
qsort(pts, NP, sizeof(point_t), cmp_x);
printf("min: %g; ", sqrt(closest(pts, NP, &a, &b)));
printf("point (%f,%f) and (%f,%f)\n", a->x, a->y, b->x, b->y);
return 0;
}</lang>
--[[User:Ledrug|Ledrug]] 21:34, 17 June 2011 (UTC)
: Actually, never mind--the wiki algorithm is bunk if implemented at face value, it fails miserably in some pathological cases. --[[User:Ledrug|Ledrug]] 22:56, 17 June 2011 (UTC)
|