/*
simplicity.c: draw "simple" numbers on the complex plane
Charlie Loyd, 2009-01-20ish.
I hereby cram this in the public domain.

Idea: http://reenigne.livejournal.com/114987.html
More examples: http://basecase.org/env/arabesques

Produces a P3-style PPM, the easiest format I could find. ImageMagick (e.g., convert) knows how to read it.

I'm aware that this is not idiomatic C. I wasn't planning on distributing it, and I just wanted it to compile and run fast and fit in one screen.

I've been using:
$ gcc-mp-4.3 -funsafe-math-optimizations -O3 -march=core2 -std=c99 simplicity-fast.c -o simp && ./simp | convert - pretty.tiff && open -a preview pretty.tiff

To take advantage of two cores, you can comment out the srandom(), compile two versions, amp one off, and have convert merge their files.
*/

#include <stdio.h>
#include <stdlib.h>
#include <complex.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
//#define count 4000000
//#define depth 100

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]);
  
  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; } } }
  
  //srandom(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)random()) % 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)*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"); } }