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

static char *points =
   "6m2o1v7{;~E|DvCq:k6m'bhbkcqhsmurnnhkdbdbh'L/U/T)a+n-q7sFtM{NtVga;M8A55@/L/'DM>H:A?>D;LCRJWPWU`X`Y`[_]h`bj^mSuFt>j3ZDPDM':r;s<t>t?q@oEoFnElDk>k<n:r'bmdmelfkgkiminhpdrbraqanbm'bmcncqbr]r[oYlZj[i]i^k_mbm''",
   *p, *image, *w, *cave, *space = "\n ";
static int *offsets, *range,
           width, height, edge,
           x, y, t, count, ax, bx, dx, dy, cave_area, width_1;
static double ft, scale;

/* Return a random number in the range of [0,1]. */
static double R()
{
   return (double)rand() / RAND_MAX;
}

/* Return scaled value. */
static double P(int index)
{
   return scale * (p[index] - 39);
}

/* 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 Interpolate(
             Interpolate(
                Interpolate(P(0), P(2)),
                Interpolate(P(2), P(4))),
             Interpolate(
                Interpolate(P(2), P(4)),
                Interpolate(P(4), P(6))));
}

/* Skip forward to the next wall character, returns 0 on success. */
static int SkipToNextWall()
{
   /* Skip forward, then check bounds. without checking bounds.
      This is guaranteed to terminate because we added a terminating
      NUL character before doing the floodfill passes.               */
   return (w += strspn(w, space)) - cave > cave_area - 4;
}

/* Generate mushroom to a particular size, returns 0 on success. */
static int GenerateMushroom()
{
   /* Fill background with black pixels.  Here we always fill the
      full buffer, and not just the part that will be drawn.  This is
      because we reduce the image size in single pixel increments, but
      image is encoded in 6 pixel rows, so there might be some leftover
      pixels from the previous iteration.                               */
   memset(image, 0, edge * edge);

   scale = dx / 90.0;
   dy = 0;

   /* Regarding p+=3: at the end of the loop, 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.                        */
   for(p = points; P(0) > 0; p += 3)
   {
      /* 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.                                */
      for(memset(offsets, 0, dx * 2 * sizeof(int)); P(2) > 0; p += 6)
      {
         for(t = 0; t < dx * 4; t++)
         {
            ft = t / (dx * 4.0);

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

            /* (x,y) are guaranteed to be within bounds due to how
               we defined our curve points.                        */
            range = offsets + y * 2;
            if( *range )
            {
               *range = *range < x ? *range : x;
               range++;
               *range = *range > x ? *range : x;
            }
            else
            {
               *range = range[1] = x;

               /* Keep track of last scanline.  We will draw only
                  up to the bottomost scanline, and drop all pixels
                  beyond this last scanline.  This saves a bit of
                  space at the end.                                 */
               dy = dy > y ? dy : y;
            }
         }
      }

      /* Draw scanlines. */
      range = offsets;
      for(y = 0; y < dx; y++, range += 2)
         for(x = *range; *range && x <= range[1]; x++)
            image[y * dx + x] = p < points + 73 ? (x + y) & 1 : 1;
   }

   /* Convert to bitmap data to sixel. */
   w = image;
   for(y = 0; y <= dy; y += 6)
   {
      /* p points at next position for writing sixels. */
      p = w;
      for(x = 0; x < dx; x++)
      {
         if( (*w++ = 0x3f +
                     image[y * dx + x] +
                     (image[(y + 1) * dx + x] << 1) +
                     (image[(y + 2) * dx + x] << 2) +
                     (image[(y + 3) * dx + x] << 3) +
                     (image[(y + 4) * dx + x] << 4) +
                     (image[(y + 5) * dx + x] << 5)) > 0x3f )
            /* Keep track of last non-blank pixel. */
            p = w;
      }

      /* Drop trailing blank pixels. */
      w = p;
      *w++ = '-';
   }
   *w = '\0';

   /* Replace wall characters with sixel data.

      Start with header:

         \x1bPq (3 bytes)
            Enter sixel mode.

            On terminals that don't recognize the escape sequence
            (e.g. windows console, putty), the escape character
            and the one following it might be dropped, and we get
            the rest of the sixel data in the output, which is enough
            to tell that it's an ASCII art cave.

            On terminals that recognize but don't support sixels
            (e.g. default xterm), we might get completely blank output.

            Only on terminals with sixel enabled (e.g. "Xterm -ti 340",
            mintty) will we get a full image.

         "1;1;%d;%d (at least 10 bytes)
            Set 1:1 aspect ratio and image size.  We need to specify
            image size ahead of time, otherwise the transparent pixels
            on the right and bottom sides are silently dropped.

            Declaring image size also causes the background to be filled.
            This is nice for terminals with white backgrounds (e.g.
            "xterm -ti 340 -rv"), where the white pixels might not be
            visible otherwise.

         #1;2;100;100;100 (16 bytes)
            Set color 1 to be RGB (100%,100%,100%).

         #1 (2 bytes)
            Set current color to color 1.

      Minimum header size is 31 bytes, but because we start test fitting
      images at larger sizes (with 3 digit edge lengths), we can expect
      an initial header size of 33 bytes.  We enforce a minimum cave
      size at the beginning of the program to ensure that we have enough
      contiguous bytes at the beginning to hold the full header.         */
   w = cave + sprintf(cave, "\x1bPq\"1;1;%d;%d#1;2;100;100;100#1", dx, dx);

   /* sprintf replaced a wall character with NUL, so here we restore
      it to what it was.                                             */
   *w = '#';
   for(p = image; *p; *w++ = *p++)
   {
      if( SkipToNextWall() )
         return 1;

      /* Check for repeated pixels, up to 99.
         Note that end of row sequences ('-') are never compressed. */
      for(dy = 0; *p - '-' && p[dy] == *p && dy < 99; dy++) {}
      if( dy > 3 )
      {
         /* Got a long enough run of repeated pixels. */
         x = strcspn(w, space);
         if( x > 3 && dy > 9 )
         {
            /* Got enough wall characters to use 4 byte sequences. */
            *w++ = '!';
            *w++ = dy / 10 + '0';
            *w++ = dy % 10 + '0';
            p += dy - 1;
         }
         else if( x > 2 )
         {
            /* Got enough wall characters to use 3 byte sequences. */
            if( dy > 9 )
               dy = 9;
            *w++ = '!';
            *w++ = '0' + dy;
            p += dy - 1;
         }
      }
   }

   /* Add padding. */
   for(; w - cave < cave_area - 3; *w++ = '$')
   {
      if( SkipToNextWall() )
         return 1;
   }

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

   /* Success. */
   return 0;
}

/* 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])) < 36 || (height = atoi(argv[2])) < 12 )
   {
      return printf("%s {width>=36} {height>=12}\n", *argv);
   }

   /* Allocate output buffer for cave data, with each row padded 1 to
      account for newline characters.

      Here we allocate a buffer that's twice the requested height.
      We will generate a cave with this double height, and then
      shrink it by half later.  This is to account for the taller
      character cells in terminals.  If we don't do this, we end up
      with caves that are more elongated vertically.  With this change,
      we will end up with more horizontally orientated caves (because
      terminal character aspect ratios tend to be just short of double
      height), but it looks generally nicer that way.                   */
   width_1 = width + 1;
   cave_area = width_1 * height;
   cave = malloc(cave_area * 2);

   /* Allocate buffer for image data, assuming that the maximum number
      of pixels to be encoded is width*height*6*24.  The 6 multiplier
      is for 6 pixels per character, and 24 multiplier is for run-length
      encoding compression factor.  Both are overly optimistic, in practice
      about half of the cave data will end up as wall characters that are
      eligible for encoding, and compression ratio is about 1:2, so a
      size of width*height would suffice for the average case.              */
   for(edge = 6; edge * edge < width * height * 6 * 24; edge += 6) {}
   image = malloc(edge * edge);

   /* Allocate buffer to store offsets.  We need edge*2 number of ints
      for the rasterization process, and width*height number of ints for
      floodfill.  width*height is larger than edge*2 for most sizes
      except the smallest size.  To account for small caves, we allocate
      twice of what's needed.                                            */
   offsets = malloc(cave_area * 2 * sizeof(int));

   if( !cave || !image || !offsets )
   {
      puts("Out of memory");
      return 1;
   }
   srand(time(NULL));
   srand(7);

   /* Fill cave with random bits. */
   p = cave;
   for(y = 0; y < height * 2; 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.

            Note that the right wall is 2 cells thick.  The extra
            cell is to account for the newline characters.        */
         *p++ = y && y < height * 2 - 1 && x && x < width - 1
            ? R() <= 0.47 ? 1 : 0
            : 15;
   }

   /* Apply smoothing function. */
   for(t = 0; t < 3; t++)
   {
      for(y = 1; y < height * 2 - 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 < 2; dy++)
            {
               for(dx = -1; dx < 2; dx++)
                  count += ReadCell();
            }
            cave[y * width_1 + x] |= count > 4 ? 2 : 0;
         }
      }

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

   /* Shrink cave to the requested size. */
   p = cave;
   for(y = 0; y < height * 2; y += 2)
   {
      for(x = 0; x < width_1; x++)
         *p++ = cave[y * width_1 + x] | cave[(y + 1) * width_1 + x];
   }

   /* Add terminating NUL character.  This is needed to make strlen work. */
   *p = '\0';

   /* Connect disjoint regions. */
   for(t = 0; t < 4; t++)
   {
      /* Find a random starting point. */
      x = R() * cave_area;
      x += strlen(cave + x);
      if( x >= cave_area )
         x = strlen(cave);
      if( x >= cave_area )
         /* Couldn't find a random starting point, which means all
            disjoint regions have been filled.                     */
         break;
      if( t )
      {
         /* If we managed to find a new starting point after the first
            pass, it means we found a new disjoint region.  Connect it
            to the previously filled region by digging a path to the
            previous starting point.

            We will stop as soon as we reach a previously filled region.
            We could remove that check and save a few bytes, and the
            end result will be some excessively long tunnels that cross
            multiple regions.                                            */
         dy = dx % width_1;
         for(; x % width_1 - dy && cave[x] - '\2'; x += x % width_1 > dy ? -1 : 1)
            cave[x] = '\0';
         for(; x - dx && cave[x] - '\2'; x += x > dx ? -width_1 : width_1)
            cave[x] = '\0';
         cave[x] = '\0';
      }

      /* Floodfill this open region.

         Floodfill is implemented by maintaining an explicit stack of
         offsets in a heap.  We could have also done this with a
         recursive function, but that has a tendency to run out of stack
         space if the requested cave size is too large.                  */
      *offsets = x;
      for(y = 1; y > 0;)
      {
         if( !cave[dy = offsets[--y]] )
         {
            /* Cave cells are encoded in three states:
               0 = open, unfilled.
               1 = wall.
               2 = open, filled.

               Here we are converting a previously open unfilled region
               to a filled region.                                      */
            cave[dy] = '\2';
            offsets[y++] = dy - width - 1;
            offsets[y++] = dy - 1;
            offsets[y++] = dy + width_1;
            offsets[y++] = dy + 1;
         }
      }

      /* Record previous starting point. */
      dx = x;
   }

   /* Convert cave data to printable characters.

      Any space not covered by a floodfill pass will be converted to walls,
      since they are disconnected from other spaces.                        */
   p = cave;
   for(y = 0; y < height; y++, *p++ = '\n')
   {
      for(x = 0; x < width; x++, p++)
         *p = *p - '\2' ? '#' : ' ';
   }

   /* Binary search to find the largest mushroom image that will
      fit cave walls.                                            */
   ax = 6;
   bx = edge;
   while( bx - ax > 1 )
   {
      dx = (ax + bx) / 2;
      if( GenerateMushroom() )
         bx = dx;  /* Failed, drop upper bound. */
      else
         ax = dx;  /* Success, raise lower bound. */
   }

   /* At this point, the upper and lower bounds differ by 1, but only
      one of them is the solution.  If we broke out of the previous loop
      by raising the lower bound, the variables are as follows:
         dx = solution (previous ax-1); ax = dx; bx = dx + 1

      If we broke out of the loop by dropping the upper bound, the
      variables as follows:
         ax = dx - 1 = solution; dx = failure (previous bx-1); bx = dx

      We will regenerate mushroom with the lower bound for the latter case. */
   if( bx == dx )
   {
      dx = ax;
      GenerateMushroom();
   }

   puts(cave);
   return 0;
}