#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <time.h>
#include <assert.h>

/* The size of the image, and a pointer to the dynamically allocated memory
 * which holds the pixel values (0 for background, 1 for crystal).  They
 * are only global to save the bother of passing them around. */
int width, height;
char *img;

/* Mark a pixel as part of the crystal. */
void
plot (int x, int y)
{
   assert (x >= 0);
   assert (x < width);
   assert (y >= 0);
   assert (y < height);

   img[x + y * width] = 1;
}

/* Return the value of the pixel at the specified position. */
int
pixel (int x, int y)
{
   assert (x >= 0);
   assert (x < width);
   assert (y >= 0);
   assert (y < height);

   return img[x + y * width];
}

/* Write the image data to the filehandle 'f', in the binary variant of the
 * PBM format - one bit per pixel. */
void
write_pbm (FILE *f)
{
   int x, y;
   int data, bits;

   /* Write the header. */
   fprintf(f, "P4\n%d %d\n", width, height);

   /* Write the binary image data. */
   for (y = 0; y < height; ++y) {
      data = bits = 0;

      for (x = 0; x < width; ++x) {
         data = (data << 1) | (pixel(x, y) ? 0 : 1);

         if (++bits == 8) {
            fputc(data, f);
            data = bits = 0;
         }
      }

      /* Write leftovers for a line if the width isn't a multiple of 8. */
      if (bits) {
         fputc(data << (8 - bits), f);
      }
   }
}

/* Return a random number in the range 0 to n-1. */
int
rand_over_range (int n)
{
   double dummy;
   return (int) (modf((double) rand() / RAND_MAX, &dummy) * n);
}

/* Set the values pointed to by 'x' and 'y' to a random location along the
 * side of the image. */
void
random_start_position (int *x, int *y)
{
   /* Decide which side to put the new particle on.  The value of 'side'
    * is 0 for the top of the image, 1 for the right, 2 for the bottom
    * or 3 for the left. */
   int side = rand_over_range(4);

   *x = (side == 0 || side == 2) ? rand_over_range(width)
                                 : (side == 1) ? width - 1 : 0;
   *y = (side == 1 || side == 3) ? rand_over_range(height)
                                 : (side == 2) ? height - 1 : 0;
}

/* Return a true value if the specified point is adjacent to any part of the
 * growing crystal. */
int
next_to_crystal (int x, int y)
{
   return (x > 0 && pixel(x - 1, y)) ||
          (x < width - 1 && pixel(x + 1, y)) ||
          (y > 0 && pixel(x, y - 1)) ||
          (y < height - 1 && pixel(x, y + 1));
}

/* Randomly change the position given by one pixel in any of the four
 * cardinal directions.  For simplicity (and because it shouldn't affect
 * the crystal's growth) we allow particles to wrap around to the other side
 * of the image. */
void
shift_point (int *x, int *y)
{
   int dir = rand_over_range(4);

   if (dir == 0)
      --*y;
   else if (dir == 1)
      ++*x;
   else if (dir == 2)
      ++*y;
   else
      --*x;

   if (*x < 0)
      *x = width - 1;
   else if (*x == width)
      *x = 0;
   else if (*y < 0)
      *y = height - 1;
   else if (*y == height)
      *y = 0;
}

/* Simulate a single particle starting at the edge of the image and moving
 * about randomly until it bumps into part of the crystal. */
void
do_particle (void)
{
   int x, y;

   random_start_position(&x, &y);

   while (!next_to_crystal(x, y))
      shift_point(&x, &y);

   plot(x, y);
}


int
main (int argc, char **argv)
{
   int iterations;
   int i;

   /* Make sure we get different random numbers each time we run this. */
   srand(time(0));

   /* Check and decode the command line arguments. */
   if (argc != 4) {
      fprintf(stderr, "Usage: %s WIDTH HEIGHT ITERATIONS > out.pbm\n",
              argv[0]);
      return 1;
   }

   width = atoi(argv[1]);
   height = atoi(argv[2]);
   iterations = atoi(argv[3]);

   /* Allocate memory to store the image in. */
   img = malloc(width * height);
   if (!img) {
      fprintf(stderr, "%s: error allocating memory for image: %s\n",
              argv[0], strerror(errno));
      return 2;
   }

   /* Clear out the image memory, setting it all to black. */
   for (i = 0; i < width * height; ++i)
      img[i] = 0;

   /* Start with a single pixel set in the middle of the image. */
   plot(width / 2, height / 2);

   /* Simulate the particles as many times as requested. */
   for (i = 0; i < iterations; ++i) {
      do_particle();

      if (i && (iterations < 100 || i % (iterations / 100) == 0))
         fprintf(stderr, "Done %d/%d particles (%d%%).\n",
                 i, iterations, (int) ((i * 100.0) / iterations + 0.5));
   }

   write_pbm(stdout);

   free(img);
   return 0;
}