#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <fcntl.h>
#include <unistd.h>
#include <pthread.h>

_Thread_local unsigned long long sum = 0;
_Thread_local unsigned long long nb_points = 0;
_Atomic unsigned long long g_sum = 0;
_Atomic unsigned long long g_nb_points = 0;
int niter = 10000000;
int nepochs = 40;
int nb_threads = 4;

void monte_carlo_pi() {
  struct drand48_data rand_buffer;
  srand48_r(pthread_self(), &rand_buffer);
  for(int i=0; i<niter; i++) {
    double x, y;
    drand48_r(&rand_buffer, &x);
    drand48_r(&rand_buffer, &y);
    double distance = x*x + y*y;
    if(distance < 1)
      sum++;
    nb_points++;
  }
}

void* thread_function(void*arg ){
  int n = nepochs / nb_threads ;
  for(int i=0; i<n; i++) {
    monte_carlo_pi();
    g_sum += sum;
    g_nb_points += nb_points;
    sum = 0;
    nb_points = 0;
  }

  return NULL;
}

int main(int argc, char**argv) {
  if(argc < 1) {
    fprintf(stderr, "Usage: %s\n", argv[0]);
    return EXIT_FAILURE;
  }

  pthread_t tids[nb_threads];
  for(int i = 0; i <nb_threads; i++) {
    pthread_create(&tids[i], NULL, thread_function, NULL);
  }

  for(int i = 0; i <nb_threads; i++) {
    pthread_join(tids[i], NULL);
  }

  double pi = (double)4*g_sum / ((double)g_nb_points);
  printf("%llu: pi = %lf\n", g_nb_points, pi);
  return EXIT_SUCCESS;
}