Closest-pair problem/C

From Rosetta Code

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
 
typedef struct {
  double x, y;
} point_t;
 
inline double distance(point_t *p1, point_t *p2)
{
  return sqrt((p1->x - p2->x)*(p1->x - p2->x) +
	      (p1->y - p2->y)*(p1->y - p2->y));
}

We need an ad hoc sort function; here is a modified version of what you can find at Quicksort.

typedef int(*comparator)(const void *, const void *);
void quick(point_t *P, int *left, int *right, comparator compar)
{
  if (right > left) {
    int pivot = left[(right-left)/2];
    int *r = right, *l = left;
    do {
      while (compar(&P[*l], &P[pivot]) < 0) l++;
      while (compar(&P[*r], &P[pivot]) > 0) r--;
      if (l <= r) {
        int t = *l;
        *l++ = *r;
        *r-- = t;
      }
    } while (l <= r);
    quick(P, left, r, compar);
    quick(P, l, right, compar);
  }
}
void m_qsort(point_t *P, int *array, int length, comparator compar)
{
  quick(P, array, array+length-1, compar);
}
 
 
int sort_index_by_x(const void *ra, const void *rb) {
  const point_t *a = ra, *b = rb;
  double d = a->x - b->x;
  return (d<0) ? -1 : ( (d==0) ? 0 : 1);
}
 
int sort_index_by_y(const void *ra, const void *rb) {
  const point_t *a = ra, *b = rb;
  double d = a->y - b->y;
  return (d<0) ? -1 : ( (d==0) ? 0 : 1 );
}

double closest_pair_simple_(point_t *P, int *pts, int num, int *pia, int *pib) {
  int i, j;
  if ( num < 2 ) return HUGE_VAL;
  double ird = distance(&P[pts[0]], &P[pts[1]]);
  *pia = pts[0]; *pib = pts[1];
  for(i=0; i < (num-1); i++) {
    for(j=i+1; j < num; j++) {
      double t;
      t = distance(&P[pts[i]], &P[pts[j]]);
      if ( t < ird ) {
	ird = t;
	*pia = pts[i]; *pib = pts[j];
      }
    }
  }
  return ird;
}

This is the entry point for the simple algorithm.

double closest_pair_simple(point_t *P, int num, int *pia, int *pib) {
  int *pts, i;
  double d;
 
  pts = malloc(sizeof(int)*num);  assert(pts != NULL);
  for(i=0; i < num; i++) pts[i] = i;
  d = closest_pair_simple_(P, pts, num, pia, pib);
  free(pts);
  return d;
}

double closest_pair_(point_t *P, int *xP, int *yP, int N, int *iA, int *iB)
{
  int i, k, j;
  int *xL, *xR, *yL, *yR, *yS;
  int lA, lB, rA, rB, midx;
  double dL, dR, dmin, xm, closest;
  int nS;
 
  if ( N <= 3 ) return closest_pair_simple_(P, xP, N, iA, iB);
 
  midx = ceil((double)N/2.0) - 1;
 
  xL = malloc(sizeof(int)*(midx+1));       assert(xL != NULL);
  xR = malloc(sizeof(int)*(N-(midx+1)));   assert(xR != NULL);
  yL = malloc(sizeof(int)*(midx+1));       assert(yL != NULL);
  yR = malloc(sizeof(int)*(N-(midx+1)));   assert(yR != NULL);
 
  for(i=0;i <= midx; i++) xL[i] = xP[i];
  for(i=midx+1; i < N; i++) xR[i-(midx+1)] = xP[i];
 
  xm = P[xP[midx]].x;
 
  for(i=0, k=0, j=0; i < N; i++) {
    if ( P[yP[i]].x <= xm ) {
      yL[k++] = yP[i];
    } else {
      yR[j++] = yP[i];
    }
  }
 
  dL = closest_pair_(P, xL, yL, midx+1, &lA, &lB);
  dR = closest_pair_(P, xR, yR, (N-(midx+1)), &rA, &rB);
 
  if ( dL < dR ) {
    dmin = dL;
    *iA = lA; *iB = lB;
  } else {
    dmin = dR;
    *iA = rA; *iB = rB;
  }
 
  yS = malloc(sizeof(int)*N);    assert(yS != NULL);
  nS = 0;
  for(i=0; i < N; i++) {
    if ( fabs(xm - P[yP[i]].x) < dmin ) {
      yS[nS++] = yP[i];
     }
  }
 
  closest = dmin;
  if (nS > 1) { 
    for(i=0; i < (nS-1); i++) {
      k = i + 1;
      while( (k < nS) && ( (P[yS[k]].y - P[yS[i]].y) < dmin ) ) {
	double d = distance(&P[yS[i]], &P[yS[k]]); 
	if (  d < closest ) {
	  closest = d;
	  *iA = yS[i];
	  *iB = yS[k];
	}
	k++;
      }
    }
  }
 
  free(xR); free(xL);
  free(yL); free(yR); free(yS);
  return closest;
}

This is the entry point for the divide&conquer algorithm.

double closest_pair(point_t *P, int N, int *iA, int *iB)
{
  int *xP, *yP, i;
  double d;
 
  xP = malloc(sizeof(int)*N); assert(xP != NULL);
  yP = malloc(sizeof(int)*N); assert(yP != NULL);
 
  for(i=0; i < N; i++) {
    xP[i] = yP[i] = i;
  }
  m_qsort(P, xP, N, sort_index_by_x);
  m_qsort(P, yP, N, sort_index_by_y);
  d = closest_pair_(P, xP, yP, N, iA, iB);
  free(xP); free(yP);
  return d;
}

Testing

#define NP 10000
 
int main()
{
  point_t *points;
  int i;
  int p[2];
  double c;
 
  srand(31415);
 
  points = malloc(sizeof(point_t)*NP);
 
  for(i=0; i < NP; i++) {
    points[i].x = 20.0*((double)rand()/(RAND_MAX+1.0)) - 10.0;
    points[i].y = 20.0*((double)rand()/(RAND_MAX+1.0)) - 10.0;
  }
 
  c = closest_pair_simple(points, NP, p, p+1);
  printf("%lf %d %d (%lf)\n", c, p[0], p[1], distance(&points[p[0]], &points[p[1]]));
  c = closest_pair(points, NP, p, p+1);
  printf("%lf %d %d (%lf)\n", c, p[0], p[1], distance(&points[p[0]], &points[p[1]]));
 
  free(points);
  return EXIT_SUCCESS;
}

The divide&conquer algorithm gave 0.01user 0.00system 0:00.11elapsed, while the brute-force one gave 1.83user 0.00system 0:01.88elapsed.