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

typedef int _;
typedef double f;

png_image I;
png_bytep z[2];
f p, q, P, Q,
  a, b, c, d, i, j, k,
  m, n,
  e, E;

_ D,o,R,O,T,H,Y[2][2], x, y, w, h[2][2];

_ g(_ u, char *v)
{
   return printf(u ? u - 1 ? "Error writing %s\n"
                           : "Error reading %s\n"
                   : "%s <input.png> <output1.png> <output2.png>\n",
                 v);
}

f r(f v)
{
   return (f)rand() / RAND_MAX * v * 22 / 7;
}

void t(f u, f v)
{
   m = a * u - b * v + c;
   n = b * u + a * v + d;
}

_ s(f u, f v, f U, f V)
{
   u -= U;
   v -= V;

   p = i * U + j * V + k;
   return u * u + v * v < p * p;
}

_ main(_ u, char **v)
{
   if( u - 4 )
      return g(0, *v);

   I.version = PNG_IMAGE_VERSION;
   if( png_image_begin_read_from_file(&I, v[1]) )
   {
      w = I.width;
      H = I.height;
      for(I.format = PNG_FORMAT_RGBA; o < 2;)
         z[o++] = (png_bytep)malloc(D = w * H * 4);
      T = *z && z[1] &&
          png_image_finish_read(&I, NULL, *z, 0, NULL);
   }
   if( !T )
      return g(1, v[1]);

   memcpy(z[1], *z, D);

   srand(time(NULL));

   p = r(29) + 20;
   p /= w > H ? w : H;

   a = cos(q = r(1)) * p;
   b = sin(q) * p;
   c = r(1);
   d = r(1);

   i = cos(q = r(2));
   j = sin(q);
   k = 0;

   for(; O < 4; O++)
   {
      t(O % 2 * w, (O & 2) / 2 * H);

      p = i * m + j * n + k;
      e = !O || e > p ? p : e;
      E = !O || E < p ? p : E;
   }

   i *= p = 0.5 * sqrt(2) / (E - e);
   j *= p;

   for(k = -e * p; x < 2; x++)
      for(y = 0; y < 2; h[x][y++] = r(2) + 2)
         Y[x][y] = r(3);

   for(y = o = 0; y < H; y++)
      for(x = 0; x < w; x++, o++)
      {
         t(x, y);
         if( sizeof('c') > 1 )
         {
            D = s(m, n, floor(m), floor(n)) ||
                s(m, n, floor(m), ceil(n)) ||
                s(m, n, ceil(m), floor(n)) ||
                s(m, n, ceil(m), ceil(n));
         }
         else
         {
            p = floor(m);
            q = floor(n);
            P = m - p;
            Q = n - q;
            if( P + Q < 0.25 ||
                P + Q > 1.75 ||
                Q - P > 0.75 ||
                P - Q > 0.75 )
            {
               P = round(m);
               Q = round(n);
               T = (_)(P + Q) & 1;

               for(O = D = 0; O < 2; D ^= R > h[O++][T] / 2)
                  if( (R = fmod(Y[O][T] + P - O * Q, h[O][T])) < 0 )
                     R += h[O][T];

               D ^= m - P > fabs(n - Q) ||
                    m - P < -fabs(n - Q);
            }
            else
               D = (_)(p + q) & 1;
         }
         memset(z[D] + o * 4, 0, 4);
      }

   for(o = 0; o < 2; o++)
      if( !png_image_write_to_file(&I, v[o + 2], 0, z[o], 0, NULL) )
         return g(2, v[o + 2]);
   return 0;
}