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

/* DEFINITIONS */
/* define TRUE and FALSE */
#define	TRUE	1
#define	FALSE	0
/* define the mathematical operations */
#define	ADD	2
#define	SUB	3
#define	MUL	4
/* define the criterion for convergence */
#define	TEST_LIMIT	2.0

/* MACROS */
/* macro to return the absolute value of a number */
#define	abs(a)		( ((a)<0) ? (-(a)):(a) )
/* macro to return the length of a vector (x, y).
   The same as the modulo of the complex number x+yi */
#define	modulo(x, y)	( sqrt( ((x)*(x))+((y)*(y)) ) )

/* DATA TYPE DEFINITIONS */
/* define a data type for a complex type. */ 
typedef	struct	_complex_number_type {
	double	Real, Imag;
} ComplexNumberType;

/* FUNCTION PROTOTYPES */
void	ComplexArithmetic(ComplexNumberType n1, ComplexNumberType n2, ComplexNumberType *r, int operation);
void	MandelbrotTransform(ComplexNumberType z, ComplexNumberType c, ComplexNumberType *r);

/* main function */
void	main(int argc, char **argv)
{
	double			StartX, StartY,
				Width, Height,
				Step, real, imag;
	int			Repetitions, Iterations, PointsInside,
				i, BelongsFlag;
	ComplexNumberType 	c, z, res;

	/* Check wheather we have the correct number of arguments. If not print usage */
	if( argc != 7 ){
		fprintf(stderr, "Usage: %s StartX StartY Width Height Step Iterations\n", argv[0]);
		exit(1);
	}

	/* Read the various parameters from the command line arguments. */
	/* Read Upper Horizontal Coordinate of Field of interest */
	StartX = atof(argv[1]);
	/* Read Upper Vertical Coordinate of Field of interest */
	StartY = atof(argv[2]);
	/* Read Width of Field of interest */
	if( (Width=atof(argv[3])) <= 0.0 ){
		fprintf(stderr, "The width of the plane must be positive.\n");
		exit(1);
	}
	/* Read Height of Field of interest */
	if( (Height=atof(argv[4])) <= 0.0 ){
		fprintf(stderr, "The height of the plane must be positive.\n");
		exit(1);
	}
	/* Read the Step or the Sampling rate */
	if( (Step=atof(argv[5])) <= 0.0 ){
		fprintf(stderr, "The grid size must be positive.\n");
		exit(1);
	}
	/* Read the number of Iterations before we can safely deduce that a point
	   belongs to the plane. */
	if( (Iterations=atoi(argv[6])) <= 0 ){
		fprintf(stderr, "The number of iterations must be positive.\n");
		exit(1);
	}

	PointsInside = 0;
	/* Go around the specified field of interest. */
	for(imag=StartY;imag<StartY+Height;imag+=Step)
		for(real=StartX;real<StartX+Width;real+=Step){
			/* Pick up a point, c */
			c.Real = real;
			c.Imag = imag;
			/* Zero z */
			z.Real = 0.0;
			z.Imag = 0.0;
			BelongsFlag = TRUE;
			/* Do a test wheather this point belongs to the Mandelbrot Set. */
			for(i=0;i<Iterations;i++){
				/* Do mandelbrot transform */
				MandelbrotTransform(z, c, &res);
				z.Real = res.Real;
				z.Imag = res.Imag;
				/* Apply criterion */
				if( (abs(z.Real)>(TEST_LIMIT)) ||
				    (abs(z.Imag)>(TEST_LIMIT)) ){
					/* If the point does not belong to the Set,
					   break. There is no need for more. */
					/* Keep the repetitions to use it later */
					BelongsFlag = FALSE;
					Repetitions = i;
					break;
				}
			}
			if( BelongsFlag ) PointsInside++;

/*			if( BelongsFlag )
				printf("BELONGS: %lf + %lfi\n", z.Real, z.Imag);
			else
				printf("DOESNOT: %lf + %lfi, %d\n", z.Real, z.Imag, Repetitions);
*/
			if( BelongsFlag ) printf("255\n");
			else printf("0\n");
		}

	fprintf(stderr, "A total of (%d x %d) %d points were tested.\n",
				(int )(Width/Step),
				(int )(Height/Step),
				(int )((Width/Step)*(Height/Step)) );
	fprintf(stderr, "%d points belong to the Mandelbrot Set.\n", PointsInside);
	fprintf(stderr, "%d points do not.\n", (int )((Width/Step)*(Height/Step)) - PointsInside);
			
}




/* Function to do the MandelBrot Transform:
		z -> z^2 + c
   Is done in two steps:
	a. Calculate z^2,
	b. add starting number, c. */
/* Parameters:
	z, c: 	Complex numbers.
	r:	Address to store the result (z^2+c).

	returns nothing.
*/
void	MandelbrotTransform(ComplexNumberType z, ComplexNumberType c, ComplexNumberType *r)
{
	ComplexNumberType	z_squ;

	/* Do z*z, put it in z_squ */
	ComplexArithmetic(z, z, &z_squ, MUL);
	
	/* Do z_squ+c, put it in the returned result. */
	ComplexArithmetic(z_squ, c, r, ADD);
}
/* End MandlebrotTransform */





/* Function implementing three Arithmetic Operations:
	Addition, Substraction, Multiplication
   on Complex Numbers */
/* Parameters:
	n1, n2: Complex numbers to operate on.
	operation: The Arithmetic operation to perform. Must be one of {ADD, SUB, MUL}.
	r:	Address to put the result in.

	returns nothing.
*/
void	ComplexArithmetic(ComplexNumberType n1, ComplexNumberType n2, ComplexNumberType *r, int operation)
{
	
	switch(operation){
		case ADD: r->Real = n1.Real + n2.Real;
			  r->Imag = n1.Imag + n2.Imag;
			  break;
		case SUB: r->Real = n1.Real - n2.Real;
			  r->Imag = n1.Imag - n2.Imag;
			  break;
		case MUL: r->Real = (n1.Real)*(n2.Real) - (n1.Imag)*(n2.Imag);
			  r->Imag = (n1.Real)*(n2.Imag) + (n1.Imag)*(n2.Real);
			  break;
		default:  break;
	}
}
/* End ComplexArithmetic */