// metropolis yöntemi ile random sayı üretir it generates random numbers by using metropolis method
// compile: gcc  -lm -fopenmp random_openmp.c
// usage: ./a.out  trial# thread#
// note: process# must be square of an integer


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


// w is gaussian weight function
double w(double x, double y)
{
   double sigma ;
// sigma is standard deviation
   sigma = 0.05 ;

   return exp(-(x *x + y*y) / (2 * sigma * sigma) ) ;

}



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


   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 ;
   }


}


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



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

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

   int total_number[tnumber];  // this array holds how many number generated in each thread, and this is shared among all threads
   

   #pragma omp parallel num_threads(tnumber) 
   {
      double * local_x_array ;
      double * local_y_array ;
      int myrank = omp_get_thread_num();
      int i ;

      local_x_array = malloc( count * sizeof(double)) ;
      local_y_array = malloc( count * sizeof(double)) ;



   // define boundries for each thread
      double x_min, x_max, y_min, y_max ;

      double interval = 2.0 / sqrt(tnumber) ; 
      x_min = -1.0 + interval * (myrank % ((int) sqrt(tnumber))) ;
      x_max = x_min + interval ;
      y_min = -1.0 + interval * (myrank / (int) sqrt(tnumber)) ;
      y_max = y_min + interval ;

      int counter = 0 ; // counter holds number of random points in local_x_array and local_y_array

      srand(myrank) ; // seed each process a different number to generate different random number


      while(counter != 1)
      {

         x = x_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to x
         y = y_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to y



         if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
         { 
            local_x_array[0] = x ;
            local_y_array[0] = y ;
            counter++ ;
         }

      }


      for(i=1 ; i < count ; i++ )
      {
         generator(&x, &y);
         if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
         { 
            local_x_array[counter] = x ;
            local_y_array[counter] = y ;
            counter++ ;
         }

      }



      total_number[myrank] = counter ;

// to wait until total_number array is populated place a barrier 

#pragma omp barrier


      if (myrank == 0)
      {
         int sum = 0 ; // sum holds sum of generated numbers by each thead 
         for(i=0 ; i < tnumber ; i++)
         {
            sum = sum + total_number[i] ;
         }
      printf("sum is %d\n", sum);
         x_array = malloc(sizeof(double) * sum);
         y_array = malloc(sizeof(double) * sum);

      }

// to prevent threads to write before thread 0 create x_array, y_array place another barrier
#pragma omp barrier

   // now each array populates global x_array, y_array by using private array local_x_array, local_y_array   
   // firstly, each thread computes displacement from initial element of global array 
      int displs = 0 ;
      
      for(i=0 ; i < myrank ; i++)
      {
         displs = displs + total_number[i];
      }


      for(i=0 ; i < counter ; i++ )
      {

         x_array[displs + i] = local_x_array[i] ;
         y_array[displs + i] = local_y_array[i] ;

      }

   } //end pragma


   return 0;

}