// calculate integral using master-workers model
// compile:  mpicc integral1_mpi.c -lm
// usage: mpirun -n process# ./a.out trial#

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <mpi.h>


// w is gaussian weight function
double w(double x, double y)
{
  return exp( - ( x*x + y*y ) );
}



void generator(double * x, double * y)
{
   double xt, yt, ratio, tmp ;
   double delta = 0.5 ;


   tmp = (double)rand() / (double)RAND_MAX  ;
// update x by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
   xt = (*x) + delta * (2 * tmp - 1) ; 

   tmp = (double)rand() / (double)RAND_MAX  ;
// update y by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
   yt = (*y) + delta * (2 * tmp - 1) ; 


// compare updated x,y values with old x,y values, accept or reject updated values as new values according to weight function
   ratio = w(xt, yt ) / w(*x, *y) ;
   tmp = (double)rand() / (double)RAND_MAX  ;

   if(ratio > tmp )
   {

      *x = xt ;
      *y = yt ;
   }


}


//f is the function to be integrated
double f(double x, double y)
{


return M_PI *  ( x*x + y*y ) ;

}


int main(int argc, char* argv[])
{
   double  x, y ;
   double * x_array ;
   double * y_array ;

   double sum ;
   double total_sum ;
   double integral ;
   int trials, trials_per_p;
   int i, j ;
   int pnumber, myrank;

   MPI_Status status;
   MPI_Init(&argc,&argv);
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   MPI_Comm_size(MPI_COMM_WORLD, &pnumber);




   x = 0.0 ;
   y = 0.0 ;


   trials = atoi(argv[1]);
   trials_per_p = trials / (pnumber - 1 );

   x_array = malloc( trials_per_p * sizeof(double)) ;
   y_array = malloc( trials_per_p * sizeof(double)) ;



   if(myrank == 0) //master process generates random numbers and sends them to workers
   {
      for(j=1 ; j < pnumber ; j++)
      {
         for(i=0 ; i < trials_per_p ; i++)
         {
            generator(&x, &y);
               x_array[i] = x ;
               y_array[i] = y ;
         }

         MPI_Send(x_array, trials_per_p , MPI_DOUBLE, j, 1234 , MPI_COMM_WORLD) ; 
         MPI_Send(y_array, trials_per_p , MPI_DOUBLE, j, 5678 , MPI_COMM_WORLD) ;    

      }
   }


   if(myrank != 0) // worker processes receive random numbers from master process and calculate sum
   {

      MPI_Recv(x_array, trials_per_p, MPI_DOUBLE, 0, 1234, MPI_COMM_WORLD, &status) ;
      MPI_Recv(y_array, trials_per_p, MPI_DOUBLE, 0, 5678, MPI_COMM_WORLD, &status) ;

      sum = 0.0 ;
      for(i=0 ; i < trials_per_p ; i++)
      {

         sum += f(x_array[i], y_array[i]) ;
      } 


   }

   MPI_Reduce(&sum, &total_sum, 1, MPI_DOUBLE , MPI_SUM, 0, MPI_COMM_WORLD) ; // all calculated sums are accumulated in master process

   
   if(myrank == 0)
   {
      integral = total_sum / ( (pnumber -1) * trials_per_p ) ;
      printf("%f %f\n", integral , M_PI);
   }


   MPI_Finalize();

   return 0;

}