/* ===========================================================
   compute Pi by integrating the function 1/SQRT(1-x^2)


   OpenMP version 1 - not using parallel for
   =========================================================== */

using namespace std;

#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

double f(double a)
{
   return( 2.0 / sqrt(1 - a*a) );
}

/* =======================
   MAIN
   ======================= */

int main(int argc, char *argv[])
{
   int    N;
   double pi;
   double w;      // Width of rectangle

   /* -----
      Check command line
      ----- */
   if ( argc != 2 )
   {
      cout << "Usage: " << argv[0] << " Num_Intervals" << endl;
      exit(1);
   }

   /* -----
      Get number of intervals and number of threads
      ----- */
   N = atoi(argv[1]);
   w = 1.0/(double) N;

   pi = 0.0;

   /* ==================================================
      Calaculate pi by adding area of rectangles
      ================================================== */

   #pragma omp parallel
   {
      int    i;
      int    num_threads;
      double x;

      num_threads = omp_get_num_threads();

      for (i = omp_get_thread_num(); i < N; i = i + num_threads)
      {
         x = (i + 0.5)*w;
         #pragma omp critical
         {
         pi = pi +  w*f(x);
         }
      }

   }

   cout << "Computed Pi = " << pi << endl << endl;

}