// Sample program to estimate PI = 3.14159... by randomly generating
// points in the positive quadrant and testing whether they are
// distance <= 1.0 from the origin.  The ratio of "hits" to tries is
// an approximation for the area of the quarter circle with radius 1
// so multiplying the ratio by 4 estimates pi.
// 
// usage: omp_picalc <num_samples> [threads]
//   num_samples: int, how many sample points to try, higher gets closer to pi
//   threads: how many threads to use, to a limit more gets better performance

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

int main(int argc, char **argv) { 
  if(argc < 2){
    printf("usage: omp_picalc <num_samples> [threads]\n");
    printf("num_samples: int, how many sample points to try, higher gets closer to pi\n");
    printf("threads: how many threads to use, to a limit more gets better performance\n");
    return -1;
  }

  if(argc > 2){
    omp_set_num_threads(atoi(argv[2]));
  }

  int npoints = atoi(argv[1]);
  int total_hits=0;

  #pragma omp parallel reduction(+: total_hits) 
  { 
    int num_threads = omp_get_num_threads(); 
    int thread_id = omp_get_thread_num();
    int points_per_thread = npoints / num_threads; 
    unsigned int seed = 123456789 * thread_id;
    int i;
    for (i = 0; i < points_per_thread; i++) { 
      double x = ((double) rand_r(&seed)) / ((double) RAND_MAX);
      double y = ((double) rand_r(&seed)) / ((double) RAND_MAX);
      if (x*x + y*y <= 1.0){
        total_hits++;
      }
    } 
  } 

  double pi_est = ((double)total_hits) / npoints * 4.0;
  printf("npoints: %8d\n",npoints);
  printf("hits:    %8d\n",total_hits);
  printf("pi_est:  %f\n",pi_est);
}