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

typedef struct _node {
   double x, y;
   struct _node *next;
} node;

typedef node *nodePointer;

typedef nodePointer **hashTable;

#define NSRCH 10

/* Create an m x m hash table of n randomly generated points */
hashTable genRandomPoints ( int n, int m )
{
   int i, j, k;
   double x, y;
   hashTable H;
   nodePointer p;

   /* First create an m x m array of list headers */
   H = (nodePointer **)malloc(m * sizeof(nodePointer *));
   for (i=0; i<m; ++i) {
      H[i] = (nodePointer *)malloc(m * sizeof(nodePointer));
      for (j=0; j<m; ++j) H[i][j] = NULL;
   }

   for (k=0; k<n; ++k) {
      /* Generate the x and y coordinates randomly in [0,1) */
      x = (double)rand() / ((double)RAND_MAX + 1);
      y = (double)rand() / ((double)RAND_MAX + 1);

      /*Locate the cell for (x,y) */
      i = (int)(x * (double)m); j = (int)(y * (double)m);

      /* Insert (x,y) in the (i,j)-th list */
      p = (node *)malloc(sizeof(node));
      p -> x = x; p -> y = y; p -> next = H[i][j];
      H[i][j] = p;
   }

   return H;
}

/* Print statistics for the hash table of n points */
void getHashStatus ( hashTable H, int n, int m )
{
   int *freq, i, j, k, max, min;
   nodePointer p;

   max = -1; min = n + 1;

   /* freq[i] stores the number of lists with exactly i entries.
      In the worst case, a list may be as long as containing all
      the input points. */
   freq = (int *)malloc((n+1) * sizeof(int));
   for (i=0; i<=n; ++i) freq[i] = 0;
   for (i=0; i<m; ++i) for (j=0; j<m; ++j) {
      /* Calculate the size of the (i,j)-th list */
      k = 0; p = H[i][j];
      while (p != NULL) { ++k; p = p -> next; }
      ++freq[k];
      if (k > max) max = k;
      if (k < min) min = k;
   }
   printf("Minimum number of points in a cell = %d\n", min);
   printf("Maximum number of points in a cell = %d\n", max);
   for (k=min; k<=max; ++k) if (freq[k] > 0)
      printf("\t%d cells with %d points\n", freq[k], k);
   free(freq);
}

/* This is the search primitive in a cell with the following parameters.
   H is the m x m hash table storing n entries.
   Search needs to be carried out in the (i,j)-th list.
   The query point is (x,y).
   min stores the current minimum distance from (x,y) to an input point
   stored in (*X,*Y). If a closer point is located in the (i,j)-th cell,
   min, *X and *Y are updated. */
double searchCell ( hashTable H, int m, int i, int j, double x, double y, double min, double *X, double *Y )
{
   nodePointer p;
   double dist;

   /* If (i,j) is not a valid cell in the m x m hash table, do nothing */
   if ((i < 0) || (i >= m) || (j < 0) || (j >= m)) return min;

   /* Otherwise, find the distance between the query point and each point in
      the current cell. If a point closer than thecurrent minimum is located,
      update min, *X and *Y accordingly. */
   p = H[i][j];
   while (p != NULL) {
      dist = sqrt((x - p->x)*(x - p->x) + (y - p->y)*(y - p->y));
      if (dist < min) { min = dist; *X = p -> x; *Y = p -> y; }
      p = p -> next;
   }

   return min;
}

/* This is the search function for finding the closest point.
   H is the m x m hash table storing n elements. The query point
   is (x,y). */
void findClosest ( hashTable H, int n, int m, double x, double y )
{
   int i, j, r, s, t;
   double min = 2, X, Y;

   /* First, locate the cell (s,t) to which the query point (x,y) belongs */
   s = (int)(x * (double)m); t = (int)(y * (double)m);

   /* Search in gradually expanding squares centered at (s,t) */
   r = 0; X = Y = -1;
   while (1) {
      /* At this point, the subsquare corresponding to r-1 is already processed.
         It is now necessary to process only the outermost layer in the search
         square for the current value of r. */
      for (i=s-r; i<=s+r; ++i) min = searchCell(H,m,i,t-r,x,y,min,&X,&Y);  /* Bottom side of the square */
      for (i=s-r; i<=s+r; ++i) min = searchCell(H,m,i,t+r,x,y,min,&X,&Y);  /* Top side of the square */
      for (j=t-r+1; j<t+r; ++j) min = searchCell(H,m,s-r,j,x,y,min,&X,&Y); /* Left side of the square */
      for (j=t-r+1; j<t+r; ++j) min = searchCell(H,m,s+r,j,x,y,min,&X,&Y); /* Right side of the squre */

      /* If a point is located, stop the initial search */
      if (min < 2) break; 

      /* Otherwise increase the size of the search square, and repeat the search */
      ++r;
   }

   printf("\n+++ Search for point closest to (%7.5lf,%7.5lf)\n", x, y);
   printf("\tInitial search gives    (%7.5lf,%7.5lf)\n", X, Y);

   /* Now starts the second phase of the search */
   i = (int)((x - min) * (double)m); /* This is the leftmost i to be considered */
   while (i < s-r) { /* The (2r+1) x (2r+1) square is already searched in the first phase */
      for (j=(y-min)*m; j<=(y+min)*m; ++j) min = searchCell(H,m,i,j,x,y,min,&X,&Y);
      ++i;
   }
   i = (int)((x + min) * (double)m); /* This is the rightmost i to be considered */
   while (i > s+r) {
      for (j=(y-min)*m; j<=(y+min)*m; ++j) min = searchCell(H,m,i,j,x,y,min,&X,&Y);
      --i;
   }
   j = (int)((y - min) * (double)m); /* This is the bottommost j to be considered */
   while (j < t-r) {
      for (i=s-r; i<=s+r; ++i) min = searchCell(H,m,i,j,x,y,min,&X,&Y);
      ++j;
   }
   j = (int)((y + min) * (double)m); /* This is the topmost j to be considered */
   while (j > t+r) {
      for (i=s-r; i<=s+r; ++i) min = searchCell(H,m,i,j,x,y,min,&X,&Y);
      --j;
   }

   printf("\tThe correct result is   (%7.5lf,%7.5lf)\n", X, Y);
}

/* Free all dynamically allocated memory*/
void windUp ( hashTable H, int m )
{
   int i, j;
   nodePointer p, q;

   for (i=0; i<m; ++i) {
      for (j=0; j<m; ++j) {
         /* Free the (i,j)-th list node by node */
         p = H[i][j];
         while (p != NULL) {
            q = p -> next;
            free(p);
            p = q;
         }
      }
      free(H[i]);
   }
   free(H);
}

int main ( int argc, char *argv[] )
{
   int n, m, i;
   hashTable H;
   double x, y;

   srand((unsigned int)time(NULL));
   if (argc >= 2) n = atoi(argv[1]); else {
      printf("n = "); scanf("%d", &n);
   }

   /* Set m to the ceiling of the square root of n */
   m = (int)sqrt(n); if (n != m * m) ++m;

   /* Populate the m x m hash table with n points */
   H = genRandomPoints(n,m);

   /* Print the statistics of the hash table */
   getHashStatus(H,n,m);

   /* Loop for search of query points */
   for (i=0; i<NSRCH; ++i) {
      /* Generate the query point randomly */
      x = (double)rand()/((double)RAND_MAX+1);
      y = (double)rand()/((double)RAND_MAX+1);

      /* Call the search function */
      findClosest(H,n,m,x,y);
   }

   /* Free memory */
   windUp(H,m);

   exit(0);
}