#include <stdio.h>
#include <math.h>
#include <omp.h>

#define image_height  4096  // Size of output image 
#define image_width   4096  // Size of output image 
#define max_count    100   // max iter for each point in image 
#define max_size     4

int image[image_height][image_width];

int main()
{
  int count, pxval;
  double xinc, yinc, xstart, ystart, yend, xend, ymid, xmid;
  double tmp, zre, zim, p_imag, p_real, mxval, mnval, t;
  char c;
  FILE *fp;
  //
  //
  t = omp_get_wtime();
  //
  // Parameters of the fractal space that are displayed 
  //
  yinc = 1e-6;
  xinc = 1e-6;
  ymid = 0.08405;
  xmid = -1.2540;
  //
  ystart = ymid - (yinc * image_height/2);
  xstart = xmid - (xinc * image_width/2);
  yend   = ymid + (yinc * image_height/2);
  xend   = xmid + (xinc * image_width/2);

#pragma omp parallel private (p_real, p_imag, zre, zim, count )
  {
    printf("I am thread #%d of %d threads..\n",omp_get_thread_num(),
	   omp_get_num_threads());
    

#pragma omp for 
    for ( int i=0; i<image_height; i++ ) {
      //
      // Determine real coordinate value for this pixel 
      //
      p_imag = ystart + i* yinc;
      for ( int j=0; j<image_width; j++ ) {
	//
	// Determine imaginary coordinate value for this pixel 
	//
	p_real = xstart + j*xinc;
	//	z0 = p_real + p_imag*I;
	//	z = z0;
	zre = p_real;
	zim = p_imag;
	count = 0;
	//
	// Iterate MAX_COUNT times or until divergence 
	//
	while ((count < max_count) && (sqrt(zre*zre+zim*zim) < max_size)) {
	  //	  z = z*z + z0;
	  tmp = zre*zre - zim*zim + p_real;
	  zim = 2*zre*zim + p_imag;
	  zre = tmp;
	  count = count + 1;
	}
	//
	// Set value in output image 
	//
	pxval = count + 2;
	image[i][j] = pxval;
      }
    }
  }
  printf("Elapsed time = %.2lf\n",omp_get_wtime()-t);
  //
  // Save image in pgm (binary) format:
  //
  mnval = 1e20; mxval=-1e20;
  for (int i=0; i<image_height; i++)
    for (int j=0; j<image_width; j++ ) {
	pxval = image[i][j];
	mnval = fmin(mnval,pxval);
	mxval = fmax(mxval,pxval);
    }
  fp = fopen("mand.pgm","w");
  fprintf(fp,"P5\n");
  fprintf(fp,"%d %d\n",image_height, image_width);
  fprintf(fp,"255\n");
  for (int i=0; i<image_height; i++)
    for (int j=0; j<image_width; j++ ) {
      c = 255*(image[i][j]-mnval)/(mxval-mnval);
      fputc(c,fp);
    }
  fclose(fp);
}