#include <stdio.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>

typedef long hint;

// Compile with gcc -DITERATION=n or uncomment the #define
// The XPM image printed to stdout will be 2^n pixels on a side
//#define ITERATION 11

#define IW (1 << ITERATION)

// Monochrome:
//#define S0 "0 0 0"
//#define S1 "16384 16384 16384"
//#define S2 "32768 32768 32768"
//#define S3 "49152 49152 49152"

// Color:
//#define S0 "5000 0 0"
//#define S1 "5000 5000 0"
//#define S2 "0 5000 0"
//#define S3 "0 0 5000"

#define S0 (60000/ITERATION), 0, 0
#define S1 (60000/ITERATION), (60000/ITERATION), 0
#define S2 0, (60000/ITERATION), 0
#define S3 0, 0, (60000/ITERATION)


int r=0, g=0, b=0;

int addcolors(int newr, int newg, int newb) {
  r += newr;
  g += newg;
  b += newb;
  //fprintf(stderr, "r+=%d g+=%d b+=%d r=%d g=%d b=%d\n",
  //  newr, newg, newb, r, g, b);
}

// iw is half the real dimensions, so we can just check x < iw
int quadrant(hint x, hint y, hint *newx, hint *newy, hint iw) {
  if(x < iw) {
    *newx = x;
    if(y < iw) {
      // 00
      // x0
      *newy = y;
      return 0;
    } else {
      // x0
      // 00
      *newy = y-iw;
      return 1;
    }
  } else {
    *newx = x - iw;
    if(y < iw) {
      // 00
      // 0x
      *newy = y;
      return 3;
    } else {
      // 0x
      // 00
      *newy = y - iw;
      return 2;
    }
  }
}

void hilbertA(hint x, hint y, hint iw);
void hilbertC(hint x, hint y, hint iw);
void hilbertD(hint x, hint y, hint iw);
void hilbertU(hint x, hint y, hint iw);

void hilbertA(hint x, hint y, hint iw) {
  iw = iw >> 1;

  hint newx, newy;
  int quad = quadrant(x,y, &newx,&newy, iw);

  switch(quad) {
    case 0: addcolors(S0); break;
    case 1: addcolors(S1); break;
    case 2: addcolors(S2); break;
    case 3: addcolors(S3); break;
  }

  if(iw == 1) {
    return;
  }
 
  switch(quad) {
    case 0: return hilbertD(newx, newy,  iw);
    case 1: return hilbertA(newx, newy,  iw);
    case 2: return hilbertA(newx, newy,  iw);
    case 3: return hilbertC(newx, newy,  iw);
  }
}

void hilbertC(hint x, hint y, hint iw) {
  iw = iw >> 1;

  hint newx, newy;
  int quad = quadrant(x,y, &newx,&newy, iw);

  switch(quad) {
    case 2: addcolors(S0); break;
    case 1: addcolors(S1); break;
    case 0: addcolors(S2); break;
    case 3: addcolors(S3); break;
  }

  if(iw == 1) {
    return;
  }
 
  switch(quad) {
    case 2: return hilbertD(newx, newy,  iw);
    case 1: return hilbertA(newx, newy,  iw);
    case 0: return hilbertA(newx, newy,  iw);
    case 3: return hilbertC(newx, newy,  iw);
  }
}

void hilbertU(hint x, hint y, hint iw) {
  iw = iw >> 1;

  hint newx, newy;
  int quad = quadrant(x,y, &newx,&newy, iw);

  switch(quad) {
    case 2: addcolors(S0); break;
    case 3: addcolors(S1); break;
    case 0: addcolors(S2); break;
    case 1: addcolors(S3); break;
  }

  if(iw == 1) {
    return;
  }
 
  switch(quad) {
    case 2: return hilbertD(newx, newy,  iw);
    case 3: return hilbertA(newx, newy,  iw);
    case 0: return hilbertA(newx, newy,  iw);
    case 1: return hilbertC(newx, newy,  iw);
  }
}


void hilbertD(hint x, hint y, hint iw) {
  iw = iw >> 1;

  hint newx, newy;
  int quad = quadrant(x,y, &newx,&newy, iw);

  switch(quad) {
    case 0: addcolors(S0); break;
    case 3: addcolors(S1); break;
    case 2: addcolors(S2); break;
    case 1: addcolors(S3); break;
  }

  if(iw == 1) {
    return;
  }
 
  switch(quad) {
    case 0: return hilbertD(newx, newy,  iw);
    case 3: return hilbertA(newx, newy,  iw);
    case 2: return hilbertA(newx, newy,  iw);
    case 1: return hilbertC(newx, newy,  iw);
  }
}


main() {
  printf("P3 %d %d 65535\n", IW, IW);
  hint x, y;
  for(x=0; x < IW; x++) {
    for(y=0; y < IW; y++) {
      r = g = b = 0;
      hilbertA(x,y, IW);
      printf("%d %d %d\n", r, g, b);
    }
  }
}