// area.c
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
struct compl {
  double real;
  double imag;
};

#define MASTER 0

int main(int argc, char *argv[]){
  int i, j, iter, numoutside = 0, numout;
  int npoints, maxiter;
  double area, error, ztemp;
  struct compl z, c;
  int rank, size, chunk, start, end;

  MPI_Init(&argc, &argv);  
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);

  if (rank == MASTER) {
    /*// Should be done by the master only
    printf("Number of points? "); scanf("%d", &npoints);
    fflush(stdout);
    printf("Maximum number of iterations? ");  scanf("%d", &maxiter);
    fflush(stdout);
    */
    npoints = atoi(argv[1]);
    maxiter = atoi(argv[2]);
    chunk = (npoints + size - 1) / size;
  }

  // Calculate & broadcast number of points for every process
  // Broadcast maximum number of iterations
  MPI_Bcast(&chunk, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(&maxiter, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
  
  start = chunk * rank;
  end = start + chunk;

  npoints = chunk * size;

/*  Outer loops run over npoints, initialise z=c
 *  Inner loop has the iteration z=z*z+c, and threshold test
 */
  for (i=start; i<end; i++) {
    for (j=0; j<npoints; j++) {
      c.real = -2.0+2.5*(double)(i)/(double)(npoints)+1.0e-7;
      c.imag = 1.125*(double)(j)/(double)(npoints)+1.0e-7;
      z=c;
      for (iter=0; iter<maxiter; iter++){
	ztemp=(z.real*z.real)-(z.imag*z.imag)+c.real;
	z.imag=z.real*z.imag*2+c.imag; 
	z.real=ztemp; 
	if ((z.real*z.real+z.imag*z.imag)>4.0e0) {
	  numoutside++; 
	  break;
	}
      }
    }
  }

  // Reduce numoutside variable
  MPI_Reduce(&numoutside, &numout, 1, MPI_INT, MPI_SUM, MASTER, MPI_COMM_WORLD);
/*
 *  Calculate area and error and output the results in the MASTER process
 */
  if (rank == MASTER) {
      area=2.0*2.5*1.125*((double)(npoints*npoints-numout)/(double)(npoints*npoints));
      error=area/(double)npoints;

      printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n",area,error);
      
  }
      MPI_Finalize();
}