#include <stdio.h>
#include <omp.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define MIN_EXP 5
#define MAX_EXP 6  // running time is about 1 minute for MAX_EXP=8

int main(int argc,char **argv)
{
  int myid, nthreads, i;
  long long int *count, n, j;
  double PI25DT = 3.141592653589793238462643;
  double x, y, mypi, mysigma, pi, sigma, error;

  omp_set_dynamic(0);
  nthreads = omp_get_max_threads();
  count = (long long int *) malloc(nthreads*sizeof(long long int));

  n = 1;
  for (i=1;i<MIN_EXP;i++) n = n*10;
  for (i=MIN_EXP;i<=MAX_EXP;i++)
    {
      n = n*10;
      #pragma omp parallel private(myid,j,x,y)
        {
          myid = omp_get_thread_num();
          srand((unsigned)time(NULL)+myid);
          count[myid] = 0;
          for (j=0;j<n;j++)
            {
              x = rand()/((double)RAND_MAX+1);
              y = rand()/((double)RAND_MAX+1);
              if ((x*x+y*y) <= 1.0) count[myid]++;
            }
        }
      pi = 0;
      sigma = 0;
      for (j=0;j<nthreads;j++)
        {
          mypi = (double)(count[j]*4)/(double)n;
          mysigma = pow((mypi-PI25DT),2);
          pi += mypi;
          sigma += mysigma;
        }
      pi = pi/(double)nthreads;
      sigma = sqrt(sigma/(double)nthreads);
      error = pi - PI25DT;
      printf("For %d trials with n = %lld:\n",nthreads,n);
      printf("pi= %.12e, error=%.12e, sigma=%.12e\n",pi,error,sigma);
    }
  return 0;
}