#define _POSIX_C_SOURCE 200809L

#define PDF
#undef PPM

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

#ifdef PDF
#include <zlib.h>
#endif

typedef struct
{
  int x0, y0, x1, y1;
} box_t;

typedef struct
{
  double rgb[3];
  box_t box;
} rect_t;

void summed_area_table(double *s, const double *y, int w, int h)
{
  for(int j = 0; j <= h; ++j)
  for(int i = 0; i <= w; ++i)
  for(int c = 0; c < 3; ++c)
    s[3*((w+1)*j+i)+c]
      = ((i==0||j==0)?0:-s[3*((w+1)*(j-1)+i-1)+c])
      + (j==0?0:s[3*((w+1)*(j-1)+i)+c])
      + (i==0?0:s[3*((w+1)*j+i-1)+c])
      + ((i==0||j==0)?0:y[3*(w*(j-1)+i-1)+c]);
}

double summed_area(const double *s, int w, int h, box_t b, int c)
{
  (void) h;
  return
      s[3*((w+1)*b.y1+b.x1)+c]
    - s[3*((w+1)*b.y1+b.x0)+c]
    - s[3*((w+1)*b.y0+b.x1)+c]
    + s[3*((w+1)*b.y0+b.x0)+c]
    ;
}

// <https://github.com/skeeto/hash-prospector#three-round-functions>
uint32_t hash(uint32_t x)
{
  x++;
  x ^= x >> 17;
  x *= 0xed5ad4bb;
  x ^= x >> 11;
  x *= 0xac4c1b51;
  x ^= x >> 15;
  x *= 0x31848bab;
  x ^= x >> 14;
  return x;
}

double fhash(uint32_t x)
{
  return hash(x) * 2.3283064365386963e-10;
}

int main(int argc, char **argv)
{
  int quality = argc > 1 ? atoi(argv[1]) : 13;
  int seed = argc > 2 ? atoi(argv[2]) : time(0);
  srand(seed);
  int candidates = 1 << quality;
  int threads = 16;
#ifdef PDF
  int G = 65536;
  char *graphics = malloc(sizeof(*graphics) * G);
  char *graphics_end = graphics + G;
  char *pdf = malloc(sizeof(*pdf) * G);
  int pdf_bytes = 0;
  char *new_pdf = malloc(sizeof(*new_pdf) * G);
  char *new_pdf_end = new_pdf + G;
  int new_pdf_bytes = 0;
#endif

  while (1)
  {
    // read ppm
    int w = 0, h = 0;
    if (2 != scanf("P6 %d %d 255", &w, &h) || w <= 0 || h <= 0)
    {
      return ! feof(stdin);
    }
    if (fgetc(stdin) != '\n')
    {
      return 1;
    }
    int n = w * h * 3;
    unsigned char *ppm = malloc(n);
    if (! ppm)
    {
      return 1;
    }
    if (1 != fread(ppm, n, 1, stdin))
    {
      free(ppm);
      return 1;
    }

    // allocate
    int r = 0;
    char *graphics_ptr = graphics;
    graphics_ptr = stpncpy(graphics_ptr, "/DeviceRGB cs\n", graphics_end - graphics_ptr);
    double *x = malloc(sizeof(*x) * n);
    double *y = malloc(sizeof(*x) * n);
//    double *xx = malloc(sizeof(*xx) * n);
    double *xy = malloc(sizeof(*xy) * n);
    double *yy = malloc(sizeof(*yy) * n);
    int N = (w + 1) * (h + 1) * 3;
    double *X = malloc(sizeof(*X) * N);
//    double *Y = malloc(sizeof(*Y) * N);
//    double *XX = malloc(sizeof(*XX) * N);
    double *XY = malloc(sizeof(*XY) * N);
    double *YY = malloc(sizeof(*YY) * N);
    if (!x || !y || /* !xx || */ !xy || !yy || !X || /* !Y || !XX || */ !XY || !YY)
    {
      return 1;
    }

    // prepare
    for (int i = 0; i < n; ++i)
    {
      double g = ppm[i] / 255.0;
      x[i] = g * g;
//      xx[i] = x[i] * x[i];
    }
    summed_area_table(X, x, w, h);
//    summed_area_table(XX, xx, w, h);
    double y0[3];
    for (int c = 0; c < 3; ++c)
    {
      y0[c] = X[((w + 1) * h + w) * 3 + c] / (w * h);
    }
    for (int i = 0; i < n; ++i)
    {
      double g = y0[i % 3];
      y[i] = g;
      xy[i] = x[i] * g;
      yy[i] = g * g;
    }
    ++r;

#ifdef PDF
    {
      int inc = snprintf(graphics_ptr, graphics_end - graphics_ptr
        , "%.2f %.2f %.2f rg %d %d %d %d re f\n"
        , sqrt(y0[0]), sqrt(y0[1]), sqrt(y0[2])
        , 0, 0, w, h
      );
      if (inc < 0)
      {
        abort();
      }
      graphics_ptr += inc;
    }

    char pageheader[128];
    snprintf
      ( pageheader, sizeof(pageheader)
      , "3 0 obj<</Type/Page/MediaBox[0 0 %d %d]/Contents 4 0 R/Parent 2 0 R/Resources<<>>>>endobj\n"
      , w, h
      );
#endif

    while (1)
    {
#ifdef PDF
      // compress graphics stream
      char *compressed_graphics = 0;
      int compressed_graphics_length = 0;
      {
        z_stream stream;
        memset(&stream, 0, sizeof(stream));
        deflateInit(&stream, Z_BEST_COMPRESSION);
        int graphics_bytes = graphics_ptr - graphics;
        uLong bound = deflateBound(&stream, graphics_bytes);
        compressed_graphics = malloc(bound);
        if (! compressed_graphics)
        {
          abort();
        }
        stream.next_in = (Bytef *) graphics;
        stream.avail_in = graphics_bytes;
        stream.next_out = (Bytef *) compressed_graphics;
        stream.avail_out = bound;
        int z = deflate(&stream, Z_FINISH);
        if (z != Z_STREAM_END)
        {
          abort();
        }
        deflateEnd(&stream);
        compressed_graphics_length = stream.total_out;
      }

      // generate pdf
      {
        char *p = new_pdf;
        char *p_end = new_pdf_end;
        p = stpncpy(p, "%PDF-1.6\n\xc2\xb7\xc2\xb7\n", p_end - p);
        int o1 = p - new_pdf;
        p = stpncpy(p, "1 0 obj<</Type/Catalog/Pages 2 0 R>>endobj\n", p_end - p);
        int o2 = p - new_pdf;
        p = stpncpy(p, "2 0 obj<</Type/Pages/Kids[3 0 R]/Count 1>>endobj\n", p_end - p);
        int o3 = p - new_pdf;
        p = stpncpy(p, pageheader, p_end - p);
        int o4 = p - new_pdf;
        int inc = snprintf(p, p_end - p, "4 0 obj<</Filter/FlateDecode/Length %d>>stream\n", compressed_graphics_length);
        if (inc < 0)
        {
          abort();
        }
        p += inc;
        if (p + compressed_graphics_length <= p_end)
        {
          memcpy(p, compressed_graphics, compressed_graphics_length);
          p += compressed_graphics_length;
        }
        p = stpncpy(p, "\nendstream endobj\n", p_end - p);
        int startxref = p - new_pdf;
        inc = snprintf
          ( p, p_end - p
          , "xref\n0 5\n0000000000 65535 f \n%010d 00000 n \n%010d 00000 n \n%010d 00000 n \n%010d 00000 n \ntrailer<</Size 5/Root 1 0 R>>\nstartxref\n%d\n%%%%EOF\n"
          , o1, o2, o3, o4, startxref
          );
        if (inc < 0)
        {
          abort();
        }
        p += inc;
        new_pdf_bytes = p - new_pdf;
        if (new_pdf_bytes > 4096)
        {
          *graphics_ptr = 0;
          fprintf(stderr, "%d\n", r - 1);
          break;
        }
        memcpy(pdf, new_pdf, new_pdf_bytes);
        pdf_bytes = new_pdf_bytes;
      }
#endif
#ifdef PPM
      if (r > 333)
      {
        break;
      }
#endif

      // find rectangle
      //summed_area_table(Y, y, w, h);
      summed_area_table(XY, xy, w, h);
      summed_area_table(YY, yy, w, h);
      double best = 1.0 / 0.0;
      box_t best_box = { 0, 0, 0, 0 };
      #pragma omp parallel for ordered schedule(static, 1)
      for (int thread = 0; thread < threads; ++thread)
      {
        double tbest = 1.0 / 0.0;
        box_t tbest_box = { 0, 0, 0, 0 };
        for (int candidate = 0; candidate < candidates / threads; ++candidate)
        {
#define R(k) fhash(((candidate * threads + thread) * 4 + k) ^ hash(seed ^ hash(r)))
          int W = exp(log(w + 1) * R(0)), H = exp(log(h + 1) * R(1)), u = (w - W) * R(2), v = (h - H) * R(3), U = u + W, V = v + H;
#undef R
          if (0 <= u && u < U && U <= w && 0 <= v && v < V && V <= h)
          {
            box_t box = { u, v, U, V };
            double a = (U - u) * (V - v);
            double score = 0;
            for (int c = 0; c < 3; ++c)
            {
              double sx = summed_area(X, w, h, box, c);
              double sxy = summed_area(XY, w, h, box, c);
              double syy = summed_area(YY, w, h, box, c);
              score += -syy + 2 * sxy - sx * sx / a;
            }
            if (score < tbest)
            {
              tbest = score;
              tbest_box = box;
            }
          }
        }

        #pragma omp ordered
        {
          if (tbest < best)
          {
            best = tbest;
            best_box = tbest_box;
          }
        }
      }

      // draw rectangle
      double draw[3], draw2[3];
      double a = (best_box.x1 - best_box.x0) * (best_box.y1 - best_box.y0);
      for (int c = 0; c < 3; ++c)
      {
        draw[c] = summed_area(X, w, h, best_box, c) / a;
        draw2[c] = draw[c] * draw[c];
      }
      for (int j = best_box.y0; j < best_box.y1; ++j)
      for (int i = best_box.x0; i < best_box.x1; ++i)
      for (int c = 0; c < 3; ++c)
      {
        int ix = (j * w + i) * 3 + c;
        y[ix] = draw[c];
        xy[ix] = x[ix] * draw[c];
        yy[ix] = draw2[c];
      }
      ++r;

#ifdef PDF
      {
        int inc = snprintf(graphics_ptr, graphics_end - graphics_ptr
          , "%.2f %.2f %.2f rg %d %d %d %d re f\n"
          , sqrt(draw[0]), sqrt(draw[1]), sqrt(draw[2])
          , best_box.x0, h - best_box.y1
          , best_box.x1 - best_box.x0, best_box.y1 - best_box.y0
        );
        if (inc < 0)
        {
          abort();
        }
        graphics_ptr += inc;
      }
#endif

    }

#ifdef PPM
    // finish
    for (int i = 0; i < n; ++i)
    {
      double g = sqrt(y[i]);
      pgm[i] = 255 * fmin(fmax(g, 0), 1);
    }

    // write ppm
    printf("P6\n%d %d\n255\n", w, h);
    if (1 != fwrite(ppm, count, 1, stdout))
    {
      free(ppm);
      return 1;
    }
#endif

    // deallocate
    free(x);
    free(y);
//    free(xx);
    free(xy);
    free(yy);
    free(X);
//    free(Y);
//    free(XX);
    free(XY);
    free(YY);
    free(ppm);

#ifdef PDF
    fwrite(pdf, pdf_bytes, 1, stdout);
    break;
#endif
  }
  return 0;
}