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

/* Target output size in bytes, not counting the final newline. */
#define SIZE 9999

static char points[] =
{
   15,70,
   11,72, 10,79, 16,84,
   20,87, 30,85, 29,79,
   28,74, 19,68, 15,70,
   0,

   59,65,
   59,68, 60,74, 65,76,
   70,78, 75,71, 71,65,
   68,61, 59,61, 59,65,
   0,

   37,8,
   46,8,  45,2,  58,4,
   71,6,  74,16, 76,31,
   77,38, 84,39, 77,47,
   64,58, 20,38, 17,26,
   14,14, 25,8,  37,8,
   0,

   29,38,
   23,33, 19,26, 24,23,
   29,20, 37,28, 43,35,
   48,41, 48,46, 57,49,
   57,50, 57,52, 56,54,
   65,57, 59,67, 55,70,
   44,78, 31,77, 23,67,
   12,51, 29,41, 29,38,
   0,

   19,75,
   20,76, 21,77, 23,77,
   24,74, 25,72, 30,72,
   31,71, 30,69, 29,68,
   23,68, 21,71, 19,75,
   0,

   59,70,
   61,70, 62,69, 63,68,
   64,68, 66,70, 66,71,
   65,73, 61,75, 59,75,
   58,74, 58,71, 59,70,
   0,

   59,70,
   60,71, 60,74, 59,75,
   54,75, 52,72, 50,69,
   51,67, 52,66, 54,66,
   55,68, 56,70, 59,70,
   0,

   0
}, *p, *image, *w, *output, *cave;
static int width, height = 6,
           *range,
           x, y, t, count, dx, dy;
static double ft, scale;

/* Interpolate between two points. */
static double Interpolate(double a, double b)
{
   return a + (b - a) * ft;
}

/* Interpolate a single component
   https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm */
static int InterpolateComponent()
{
   return (int)Interpolate(
                  Interpolate(
                     Interpolate(scale * *p,   scale * p[2]),
                     Interpolate(scale * p[2], scale * p[4])),
                  Interpolate(
                     Interpolate(scale * p[2], scale * p[4]),
                     Interpolate(scale * p[4], scale * p[6])));
}

/* Generate mushroom to a particular size. */
static void GenerateMushroom(int size)
{
   /* Set width so that (width*height/6 + height) is roughly equal to target
      size, but not exceeding target size.

      Also try to maintain square canvas.                                    */
   for(; (width = (size - 33 - height) * 6 / height) > height; height += 6) {}

   /* Compute padding size. */
   x = size - 33 - (width + 1) * height / 6;

   /* Allocate and zero-initialize output buffer.  We need enough bytes
      for header and padding, plus one byte per pixel.  We don't need
      extra bytes for footer because the image will be converted from
      1 byte per pixel to 1/6 byte per pixel before adding the footer.  */
   output = calloc(31 + x + width * height, 1);
   range = malloc(height * 2 * sizeof(int));
   if( output == NULL || range == NULL )
   {
      puts("Out of memory");
      exit(EXIT_FAILURE);
   }
   scale = width / 90.0;

   /* Initialize output header with padding. */
   strcpy(output, "\x1bPq#0;2;0;0;0#1;2;100;100;100#1");
   for(image = output + 31; image < output + 31 + x;)
      *image++ = '$';

   p = points;
   while( *p != 0 )
   {
      /* Collect horizontal ranges for each scanline.

         This allows us to draw filled shapes one scanline at a time,
         and the shapes can be both convex and concave as long as each
         scanline is continuous.  We can enforce the continuous requirement
         by tweaking our shapes, saving us the complexity of having to
         write a more featureful rasterizer.                                */
      memset(range, 0, height * 2 * sizeof(int));
      for(; p[2] != 0; p += 6)
      {
         for(t = 0; t < height * 4; t++)
         {
            ft = (double)t / (height * 4.0);

            x = InterpolateComponent();
            p++;
            y = InterpolateComponent();
            p--;

            if( x > 0 && x < width && y >= 0 && y < height )
            {
               if( range[y * 2] == 0 )
               {
                  range[y * 2] = range[y * 2 + 1] = x;
               }
               else
               {
                  range[y * 2] = range[y * 2] < x ? range[y * 2] : x;
                  range[y * 2 + 1] = range[y * 2 + 1] > x ? range[y * 2 + 1] : x;
               }
            }
         }
      }

      /* Draw scanlines. */
      for(y = 0; y < height; y++)
      {
         if( range[y * 2] != 0 )
         {
            for(x = range[y * 2]; x <= range[y * 2 + 1]; x++)
               image[y * width + x] = p < points + 73 ? (x + y) & 1 : 1;
         }
      }

      /* p[0],p[1] is the last pair of coordinates.  p[2] is a trailing zero.
         Thus, to get to the start of the next shape, we need to skip 3.     */
      p += 3;
   }

   /* Convert to Sixel.  Each set of 6 rows costs width+1 bytes. */
   w = image;
   for(y = 0; y < height; y += 6)
   {
      for(x = 0; x < width; x++)
      {
         *w++ = 0x3f +
                image[y * width + x] +
                (image[(y + 1) * width + x] << 1) +
                (image[(y + 2) * width + x] << 2) +
                (image[(y + 3) * width + x] << 3) +
                (image[(y + 4) * width + x] << 4) +
                (image[(y + 5) * width + x] << 5);
      }
      *w++ = '-';
   }

   /* Write footer. */
   strcpy(w, "\x1b\\");
}

/* Read a single cell with offset. */
static int ReadCell()
{
   return cave[(y + dy) * (width + 1) + x + dx] & 1;
}

int main(int argc, char **argv)
{
   if( argc != 3 ||
       (width = atoi(argv[1])) < 32 || (height = atoi(argv[2])) < 16 )
   {
      return printf("%s {width} {height}\n", *argv);
   }

   /* Allocate output buffer, with each row padded 1 to account for
      newline characters.                                           */
   cave = malloc((width + 1) * height + 1);
   if( cave == NULL )
   {
      puts("Out of memory");
      return 1;
   }
   srand(time(NULL));

   /* Fill cave with random bits. */
   p = cave;
   for(y = 0; y < height; y++)
   {
      for(x = 0; x < width + 1; x++)
      {
         /* We use 1 to indicate walls and 0 to indicate spaces.
            Only the lowest bit is used in the final output, but
            we fill 4 bits for the edge cells since the smoothing
            function will shift all cells right 3 bits.           */
         *p++ = (y == 0 || y == height - 1 || x == 0 || x == width - 1)
            ? 15
            : (double)rand() / RAND_MAX <= 0.47 ? 1 : 0;
      }
   }

   /* Apply smoothing function. */
   for(t = 0; t < 3; t++)
   {
      for(y = 1; y < height - 1; y++)
      {
         for(x = 1; x < width - 1; x++)
         {
            /* Sum the cave cell values from bit 0, then apply the sum
               to bit 1.

               We could also apply the sum directly to bit 0, which will
               result in an overly smoothed cave that is biased toward
               upper left.                                               */
            count = 0;
            for(dy = -1; dy <= 1; dy++)
            {
               for(dx = -1; dx <= 1; dx++)
                  count += ReadCell();
            }
            cave[y * (width + 1) + x] |= count >= 5 ? 2 : 0;
         }
      }

      /* Shift all cells right. */
      p = cave;
      for(x = 0; x < (width + 1) * height; x++)
         *p++ >>= 1;
   }

   /* Convert cells to printable characters. */
   p = cave;
   for(y = 0; y < height; y++)
   {
      for(x = 0; x < width; x++, p++)
         *p = *p ? '#' : ' ';
      *p++ = '\n';
   }
   *p = '\0';

   puts(cave);
   return 0;
}