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

int main(int argc, char **argv)
 {
  int i;
  int n=10000;
  double PI25DT = 3.141592653589793238462643;
  double pi, h, area, x, start, end;
  if (argc < 2) 
   {
    printf("Usage omp_pi.exe <number_of_threads>\n");
    exit(0);
   }
  if (!omp_get_dynamic())
   {
    //Set the number of threads
    omp_set_num_threads(atoi(argv[1]));
   }
  else
   {
    printf("ERROR: you must set OMP_DYNAMIC=FALSE\n");
    exit(0);
   }
  /* Compute an approximation to integral[4 d(arctan(x))] from x=0 to x=1 */
  start = omp_get_wtime();
  h = 1.0 / (double) n;
  area = 0.0;
  #pragma omp parallel for private(x) reduction(+:area)
  for (i = 1; i <= n; i++)
    {
      x = h * ((double)i - 0.5);
      // #pragma omp critical not needed here due to reduction clause above
      area += (4.0 / (1.0 + x*x));
    } // for
  pi = h * area;
  end = omp_get_wtime();
  printf("pi is approximately %.16f, Error is %.16f\n",pi,fabs(pi - PI25DT));
  printf("Ran in %0.5f seconds with %d threads\n",end-start, \
         omp_get_max_threads());
  //printf("%d\t%0.5f\n",omp_get_max_threads(),end-start);
  return 0;
} // main