/*
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 400.0 // was 1000
#define max 65535
#define count 20000000
#define depth 50

int main() {
	double maxseen = 0.0;
	double scalefactor;
	complex double pt;
	double px[2];
	double new;
	double*** pixels;
	
	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; }
			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"); } }