/*

gcc -funsafe-math-optimizations -O3 -std=gnu99 Klee.c -o simp -lm

date; for i in $(seq 0 7); do ./simp 16000 200000 3$i | convert -gamma 3 - slice4/s$i.png & done

*/

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
//#include <sys/types.h>
//#include <unistd.h>
#include <math.h>

#define MIN(a, b)  (((a) < (b)) ? (a) : (b))

#define W 4000
#define H 4000

#define zoom 800.0
#define max 65535

int main(int argc, char *argv[]) {
  double maxseen = 0.0;
  double scalefactor;
  complex double pt;
  double px[2];
  double new;
  double*** pixels;
  
  int count = atoi(argv[1]);
  int depth = atoi(argv[2]);

  unsigned short seed[3];
  seed[0] = atoi(argv[3]);
  seed[1] = atoi(argv[3]);
  seed[2] = atoi(argv[3]);
  seed[3] = atoi(argv[3]);
  seed48(seed);

  pixels = malloc(W * H * 3 * sizeof(double));
  for (int y = 0; y < H; y++) {
    pixels[y] = malloc(W * 3 * sizeof(double));
    for (int x = 0; x < W; x++) {
      pixels[y][x] = malloc(3 * sizeof(double));
      for (int chan = 0; chan < 3; chan++) {
	pixels[y][x][chan] = 0.0; } } }
  
  for (int c = 0; c < count; c++) {
    int r;
    unsigned fx, fy;
    pt = 0.0 + 0.0*I;
    for (int d = 0; d < depth; d++) {
      r = ((int)nrand48(seed)) % 3;
      if (r == 0) {
	pt = pt - 1 + I; }
      else if (r == 1) {
	pt = I*pt; }
      else {
	pt = cpow(pt, I); }
      
      fy = (unsigned) lround((cimag(pt)-2)*zoom+H/2);
      if (fy % H != fy) { continue; }
      
      fx = (unsigned) lround(creal(pt)*zoom+W/2);
      if (fx % W != fx) { continue; }
      
      new = pixels[fy][fx][r];
      
      new = pow(new, 0.999) + (float)count/(float)depth/2;
      
      pixels[fy][fx][r] = new;
      if (new > maxseen) {
	maxseen = new; } } }
  
  scalefactor = maxseen / ((double) max);
  
  printf("P3\n%i %i\n%i\n", W, H, max);
  for (int y = 0; y < H; y++) {
    for (int x = 0; x < W; x++) {
      for (int chan = 0; chan < 3; chan++) {
	printf("%i ",
	       (unsigned) MIN((pixels[y][x][chan]/maxseen)*max, max)); } }
    printf("\n"); } }