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

static long iterations = 1000000;

int main(int argc, char *argv[])
{
  double x, y, z; // x,y value for the random coordinate. z is used to check if x^2+y^2<=1
  long count = 0; // count holds all the number of how many good coordinates
  double pi;      // holds approx value of pi
  int numthreads = 3;
  long i;
  double stime;

  stime = omp_get_wtime();
  unsigned short xi[3];
  xi[0] = 1;
  xi[1] = 1;
  xi[2] = stime;
#pragma omp parallel
  {
#pragma omp single
    printf(" %d threads ", omp_get_num_threads());

    // srand48((int)time(NULL) ^ omp_get_thread_num()); // Give random() a seed value
#pragma omp for reduction(+ : count) private(x, y, z)
    for (i = 0; i < iterations; ++i)
    {
      x = erand48(xi);
      y = erand48(xi);
      z = ((x * x) + (y * y));

      if (z <= 1)
        ++count;
    }
  }

  pi = ((double)count / (double)(iterations)) * 4.0;
  printf("Estimated Pi: %f\n", pi);
  printf("Time: %f\n", omp_get_wtime() - stime);
  return 0;
}
