// Program 3.20: closest points

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

#include <iostream>

float randFloat () {
  return (1.0 * rand () / RAND_MAX) ;
}

struct point {
  float x; float y;
} ;

float distance (point p1, point p2) {
  float dx = p1.x - p2.x;
  float dy = p1.y - p2.y;

  return (sqrt (dx * dx + dy * dy));
}

struct node {
  point p; node *next;
  node(point pt, node* t) {
    p = pt; next = t;
  }
} ;

typedef node *mylink;

static mylink **grid;

static int cnt = 0;
static float d;
static int G;

void
gridinsert (float x, float y) {
  int i, j;

  int X = (int) (x*G) + 1;
  int Y = (int) (y*G) + 1;

  point p; p.x = x; p.y = y;

  mylink s;
  mylink t = new node(p, grid[X][Y]);

  for (i = X-1; i <= X+1; i++)
    for (j = Y-1; j <= Y+1; j++)
      for (s = grid[i][j]; s != 0; s = s->next)
	if (distance(s->p, t->p) < d) cnt++;
  grid[X][Y] = t;
}

mylink **malloc2d (int r, int c) {
  mylink **t = new mylink*[r];

  for (int i = 0; i < r; i++)
    t[i] = new mylink[c];

  return t;
} 

#define cout std::cout
#define endl std::endl

int main (int agrc, char *argv[]) {
  int i;

  int N = atoi(argv[1]);

  d = (float) atof (argv[2]);

  G = (int) (1.0 / d);

  grid = malloc2d (G+2, G+2);

  for (int i = 0; i < G+2; i++)
    for (int j = 0; j < G+2; j++)
      grid[i][j] = 0;

  for (i = 0; i < N; i++)
    gridinsert (randFloat (), randFloat ());

  cout << cnt << " pairs within " << d << endl;

  cout << "PI = " << 2.0 * cnt / (1.0 * N * (N - 1) * d * d) << endl;
}