/*

Make fractal(?) point diffusion figures.

By Charlie Loyd, January 2010. Public domain.

Far as I know, this idea was first developed by Andrew Jenner in http://www.reenigne.org/blog/a-little-mathematical-game/ , http://www.reenigne.org/blog/simplicity-density-function/ , etc.

Some samples from this code: http://basecase.org/env/arabesques .

Improvements:

- Parallelize internally. Just fork() or whatever.
- Use proper permutations of the rules instead of randomness.
- Don't be in C.

*/

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

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

#define W 4000 // width and height of the output image
#define H 4000

// arguable bug: this is absolute, not relative to W and H
#define zoom 800.0

// number of levels in the image; 2^16 is 16-bit:
#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]); // bounds-checking? what's that?
  int depth = atoi(argv[2]);

  // 48 bits of random, done ugly-ass
  // we're using the 48 functions to avoid global random state
  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)); // C sucks
  for (int y = 0; y < H; y++) { // C sucks
    pixels[y] = malloc(W * 3 * sizeof(double)); // C sucks
    for (int x = 0; x < W; x++) { // C sucks
      pixels[y][x] = malloc(3 * sizeof(double)); // C sucks
      for (int chan = 0; chan < 3; chan++) { // C sucks
  pixels[y][x][chan] = 0.0; } } } // C sucks

  for (int c = 0; c < count; c++) {
    int r;
    unsigned fx, fy;
    pt = 0.0 + 0.0*I;

    // begin interesting part
    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); }
        // end interesting part

      fy = (unsigned) lround((cimag(pt)/*-2*/)*zoom+H/2); // the -2i region is pretty
      if (fy % H != fy) { continue; } // ignore runaways (bad idea?)
      fx = (unsigned) lround(creal(pt)*zoom+W/2);
      if (fx % W != fx) { continue; }
      
      new = pixels[fy][fx][r];
      
      // doing some internal gamma, which is definitely a terrible idea:
      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++) { // C sucks
    for (int x = 0; x < W; x++) { // C sucks
      for (int chan = 0; chan < 3; chan++) { // C sucks
  printf("%i ",
         (unsigned) MIN((pixels[y][x][chan]/maxseen)*max, max)); } }  // C sucks
    printf("\n"); } }