// calculate integral using master-workers model
// compile: gcc   -fopenmp integral1_openmp.c -lm
// usage: ./a.out  trial# thread#



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


// w is 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 = 0.0 ;

   int  count = atoi(argv[1]); // how many random number will generated 

   int tnumber = atoi(argv[2]); // tnumber is total thread number, myrank is rank of each thread

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

   int ready_thread[tnumber] ; //this array holds which thread is ready to run

   #pragma omp parallel num_threads(tnumber) reduction ( +:sum ) 
   {
      int myrank = omp_get_thread_num();
      int i, j ;
      int number_for_each = count / (tnumber-1) ; // how many random number is used in each thread, master thread is not counted

      
      if(myrank == 0) //if thread no=0 then generate random numbers, else use generated random numbers
      {

         x = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1
         y = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1      

         for(j=1 ; j < tnumber ; j++)
         {
            for(i=0 ; i < count ; i++ )
            {
               x_array[i] = x ;
               y_array[i] = y ;
               generator(&x, &y);
            }
            ready_thread[ j ] = 1;

         }
      }
      else
      {
         int start_index = (myrank - 1) * number_for_each;

         while(ready_thread[myrank] != 1)  
         {
            //hold thread until ready bit becomes 1
         }
         
         for(i=0 ; i < number_for_each ; i++ )
         {
            sum += f( x_array[start_index+i] , y_array[start_index+i] ) ;
         }

      }
      

      #pragma omp barrier
   } //end pragma

   printf("%f\n", sum / count );

   return 0;
}