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

#include "mpi.h"

int main( int argc, char **argv )

{
   int n;
   int nt;
   int ite;
   float *image; 
   float *image_g;
   double s;
   int x;
   int y;
   int i;
   int j;
   double zr;
   double zi;
   double  cr;
   double  ci;
   double tr;
   double ti;
   int flag;
   double sq_mod;
   int nproc; 
   int iproc;
   int y_first;
   int y_last;

   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &nproc);
   MPI_Comm_rank (MPI_COMM_WORLD, &iproc);
   
   if (argc != 3)
      perror("usage: mand <number of iterations> <number of grid points>");
   ite = atoi(argv[1]);
   nt = atoi(argv[2]);

   if (nt % nproc) 
      perror ("number of grid points should be divisible by number of processors");

   n = nt / nproc;

   image = (float *) malloc (nt * n * sizeof(float));
   if (!image)
     perror ("malloc image");

   image_g = (float *) malloc (nt * nt * sizeof(float));
   if (!image_g)
      perror ("malloc image_g");

   s = 4.0 / ((double) (nt - 1)); 

   i = 0 ;

   y_first = n * iproc;

   y_last = y_first + n;

   for (y = y_first ; y < y_last ; y++) 
   {

      ci = -2.0 + ((double) y) * s;

      for (x = 0 ; x < nt ; x++) 
      {
         zr = 0.0;
         zi = 0.0;
         cr = -2.0 + ((double) x) * s;

         flag = 0;

         for (j = 0 ; j < ite ; j++) 
         {
           tr = (zr * zr - zi * zi) + cr;
           ti = 2.0 * zr * zi + ci;
           zr = tr ; 
           zi = ti;

           sq_mod = zr * zr + zi * zi;

           if (sq_mod > 2.0) { flag=1; break; }

         }

         if (flag) 
         {
         image[i] = ((float) j); 
         }
         else image[i] = 0.0;
         i ++;
      }
   } 

   MPI_Gather (image, nt * n, MPI_FLOAT, image_g, nt * n, MPI_FLOAT, 
      0, MPI_COMM_WORLD);
       
   if (iproc == 0) 
   {
      FILE *fp;
      int c;
      fp = fopen ("mand.out", "w");
      if (!fp) 
         perror("Error opening file mand.out");

      c = fwrite(image_g, sizeof(float), nt * nt, fp);

      if (c != (nt * nt))
         perror ("Error writing to file mand.out");
   }

   free(image);
   free(image_g);

   MPI_Finalize();

   return 0;
}
