Talk:Closest-pair problem/C

your code does NOT RUN when i compile with devc or Turbo C???? (unsigned comment added by 113.22.126.190 at 21:06, 24 October 2010)

Correct it? Try setting compiler flags for compatibility with a specific C standard? --Michael Mol 13:11, 25 October 2010 (UTC)

"The code does not run"... is like pretending to derive solar physics by the statement "the sun shines"... By the way "devc" (Dev C++ IDE?) usually is used with (an old version of) gcc. On my test machine (GNU/Linux) it runs, except for some "evil dataset" that someone gave me once upon a time... I've inspected the code with valgrind, debugged it, ... but I was not able to unwind the flow that gives the problem... I know this code hides some oddity somewhere. Likely the better thing is to rewrite it from scratch, but I've not the courage yet! :) — ShinTakezou 17:52, 30 May 2011 (UTC)

Propose replacing code

I suggest replace current code sample with rewritten code below. Reasons: 1. It doesn't segfault with 200,000 points or more; 2. It's shorter and quite a bit faster; 3. It's cleaner IMO. <lang C>#include <stdio.h>

1. include <stdlib.h>
2. include <values.h>
3. include <math.h>
4. include <string.h>

typedef struct { double x, y; } point_t, *point;

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;

}

int cmp_x(const void *a, const void *b) {

       return cmp_dbl( (*((point*)a))->x, (*((point*)b))->x );

}

int cmp_y(const void *a, const void *b) {

       return cmp_dbl( (*((point*)a))->y, (*((point*)b))->y );

}

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;

}

double closest(point* sx, int nx, point* sy, int ny, point *a, point *b) {

       int left, right, i;
       double d, min_d, x0, x1, mid, x;
       point a1, b1;
       point *s_yy;

       if (nx <= 8) return brute_force(sx, nx, a, b);

       s_yy  = malloc(sizeof(point) * ny);
       mid = sx[nx/2]->x;

       /* adding points to the y-sorted list; if a point's x is less than mid,
          add to the begining; if more, add to the end backwards, hence the
          need to reverse it */
       left = -1; right = ny;
       for (i = 0; i < ny; i++)
               if (sy[i]->x < mid) s_yy[++left] = sy[i];
               else                s_yy[--right]= sy[i];

       /* reverse the higher part of the list */
       for (i = ny - 1; right < i; right ++, i--) {
               a1 = s_yy[right]; s_yy[right] = s_yy[i]; s_yy[i] = a1;
       }

       min_d = closest(sx, nx/2, s_yy, left + 1, a, b);
       d = closest(sx + nx/2, nx - nx/2, s_yy + left + 1, ny - left - 1, &a1, &b1);

       if (d < min_d) { min_d = d; *a = a1; *b = b1; }
       d = sqrt(min_d);

       /* get all the points within distance d of the center line */
       left = -1; right = ny;
       for (i = 0; i < ny; i++) {
               x = sy[i]->x - mid;
               if (x <= -d || x >= d) continue;

               if (x < 0) s_yy[++left]  = sy[i];
               else       s_yy[--right] = sy[i];
       }

       /* compare each left point to right point */
       while (left >= 0) {
               x0 = s_yy[left]->y + d;

               while (right < ny && s_yy[right]->y > x0) right ++;
               if (right >= ny) break;

               x1 = s_yy[left]->y - d;
               for (i = right; i < ny && s_yy[i]->y > x1; i++)
                       if ((x = dist(s_yy[left], s_yy[i])) < min_d) {
                               min_d = x;
                               d = sqrt(min_d);
                               *a = s_yy[left];
                               *b = s_yy[i];
                       }

               left --;
       }

       free(s_yy);
       return min_d;

}

1. define NP 1000000

int main() {

       int i;
       point a, b;

       point pts  = malloc(sizeof(point_t) * NP);
       point* s_x = malloc(sizeof(point) * NP);
       point* s_y = malloc(sizeof(point) * NP);

       for(i = 0; i < NP; i++) {
               s_x[i] = pts + i;
               pts[i].x = 100 * (double) rand()/RAND_MAX;
               pts[i].y = 100 * (double) rand()/RAND_MAX;
       }

/* printf("brute force: %g, ", sqrt(brute_force(s_x, NP, &a, &b)));

       printf("between (%f,%f) and (%f,%f)\n", a->x, a->y, a->x, a->y);        */

       memcpy(s_y, s_x, sizeof(point) * NP);
       qsort(s_x, NP, sizeof(point), cmp_x);
       qsort(s_y, NP, sizeof(point), cmp_y);

       printf("min: %g; ", sqrt(closest(s_x, NP, s_y, NP, &a, &b)));
       printf("point (%f,%f) and (%f,%f)\n", a->x, a->y, b->x, b->y);

       /* not freeing the memory, let OS deal with it.  Habit. */
       return 0;

}</lang> --Ledrug 03:34, 18 June 2011 (UTC)