// calculate integral using client-server model
// compile: gcc   -fopenmp integral2_openmp.c -lm
// usage: ./a.out  trial# thread#
// note: trial# must be multiple of CHUNK_SIZE


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

#define CHUNK_SIZE 1000

// 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 generated = 0 ;
   int consumed = 0 ;

   #pragma omp parallel num_threads(tnumber)  
   {
      int myrank = omp_get_thread_num();
      int i, j ;
      int start_index;
      double local_sum = 0 ;
      
      if(myrank == 0 ) //if master thread 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(generated=0 ; generated < count ; generated++)
         {
            x_array[generated] = x ;
            y_array[generated] = y ;
            generator(&x, &y);

         }
      }
      else 
      {


         while( 1 )
         {

            # pragma omp critical
            {
            if( (generated - consumed) >= CHUNK_SIZE )
            {

                start_index = consumed ;
                consumed = consumed + CHUNK_SIZE ;
                for( i=0 ; i<CHUNK_SIZE ; i++)
                {
                   local_sum = local_sum + f(x_array[start_index+i], y_array[start_index+i]) ;
                }
                
                sum = sum + local_sum ;

                local_sum = 0;
            }
            } // end pragma omp critical

            if ( (generated == consumed) && (generated > 0) )
            {
               break ;
            }
          
         }

      }
      

      #pragma omp barrier
   } //end pragma

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


   return 0;
}