#include <stdio.h>

/* These determine the size of the (square) image, and total computation */
#define MAX_LOOPS 256
#define IMAGE_SIZE 256
#define IMAGE_SCALE (1.0f / (float)IMAGE_SIZE)

/* Subset of the complex plane to compute */
#define WIDTH 3.0f
#define TOP 1.5f
#define LEFT 2.0f

/* big image sizes can be too big for the stack, so must be global or on heap */
unsigned char image[IMAGE_SIZE][IMAGE_SIZE];

/* complex data type and operations */
/* put into functions to make profiling easier */
typedef struct {
    float re, im;
} Complex;

Complex Mul(Complex a, Complex b) {
    Complex result = { a.re * a.re - b.im * b.im, a.re * b.im + a.im * b.re};
    return result;
}

Complex Add(Complex a, Complex b) {
    Complex result = { a.re + b.re, a.im + b.im};
    return result;
}

float SquaredLength(Complex a) {
    return a.re * a.re + a.im * a.im;
}

/****************************************/
int main() {
    /* offset to top-left corner of image, inset half a pixel to the pixel center */
    Complex corner = { -LEFT + 0.5f * IMAGE_SCALE, -TOP + 0.5f * IMAGE_SCALE };
    FILE *imageFile;

    /* loop over image pixels */
    for(int y = 0; y < IMAGE_SIZE; ++y) {
        for(int x = 0; x < IMAGE_SIZE; ++x) {
	    /* c and initial z for this pixel */
	    Complex z = {0.0f, 0.0f};
	    Complex c = { WIDTH * (float)x * IMAGE_SCALE, WIDTH * (float)y * IMAGE_SCALE};
	    c = Add(c, corner);

	    /* iterate until too many iterations or outside 2-unit circle */
	    int i;
	    for(i=0; SquaredLength(z) < 4.0f && i < MAX_LOOPS; ++i)
		z = Add(Mul(z,z), c);

	    /* map to color for image pixel */
	    image[y][x] = 255 * (i & 0xf) / 0xf;
	}
    }

    /* pgm is an extremely simple image format to write */
    /* something like ImageMagick convert can turn it into png or tif */
    imageFile = fopen("mandel.pgm","wb");
    fprintf(imageFile, "P5\n%d %d\n255\n", IMAGE_SIZE, IMAGE_SIZE);
    fwrite(image, IMAGE_SIZE*IMAGE_SIZE, 1, imageFile);
    fclose(imageFile);

    return 0;
}