// pi_omp.c integral x=0 to 1 4.0/(1+x^2) dx
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>

#define Pi 3.1415926535897932384626433832795028841971

int main(int argc, char * argv[])
{
  long int n = 10;
  double h, pi, start;
  long int i, k;
  int numprocs, numthread = 0;

  printf("pi_omp.c running \n");
  printf("pi = integral x=0 to 1 4.0/(1+x^2) dx \n");
  printf("pi += sum h*4.0/(1.0+x*x) x=0..1  n sums\n");
  start = (double)time(NULL);
  numprocs = omp_get_num_procs();
  numthread = omp_get_thread_limit();
  printf("numprocs=%d, threadlimit=%d \n", numprocs, numthread);
  fflush(stdout);

  k = 0;
  while(k<16) // approximate digits precision (-5)
  {
    pi = 0;
    h = 1.0/(double)(n-1);
    #pragma omp parallel
    {
      double x; // local in each thread, could use  private(x)

        #pragma omp for reduction(+:pi)
        for(i=0; i<n; i++)
        {
          x = h*(double)i;
          // #pragma omp atomic  // much slower than "reduction"
          pi += h*4.0/(1.0+x*x); 
        }
    }
    #pragma omp barrier  // be sure all threads finished
    printf("n=%ld, h=%e, error=%e, took %f sec \n",
            n, h, pi-Pi, (double)time(NULL)-start);
    fflush(stdout);
    n *= 4;
    k++;
  } // end while
  printf("pi_omp.c finished \n");
  return 0;
} // end pi_omp.c