/*
  Purpose:

    This code estimates the value of PI, using Open MP.

  Discussion:

    The calculation is based on the formula for the indefinite integral:

      Integral 1 / ( 1 + X**2 ) dx = Arctan ( X ) 

    Hence, the definite integral

      Integral ( 0 <= X <= 1 ) 1 / ( 1 + X**2 ) dx 
      = Arctan ( 1 ) - Arctan ( 0 )
      = PI / 4.

    A standard way to approximate an integral uses the midpoint rule.
    If we create N equally spaced intervals of width 1/N, then the
    midpoint of the I-th interval is 

      X(I) = (2*I-1)/(2*N).  

    The approximation for the integral is then:

      Sum ( 1 <= I <= N ) (1/N) * 1 / ( 1 + X(I)**2 )

    In order to compute PI, we multiply this by 4; we also can pull out
    the factor of 1/N, so that the formula you see in the program looks like:

      ( 4 / N ) * Sum ( 1 <= I <= N ) 1 / ( 1 + X(I)**2 )

    Until roundoff becomes an issue, greater accuracy can be achieved by 
    increasing the value of N.

  Licensing:

    This code is distributed under the GNU LGPL license. 
*/
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

int main (int argc, char *argv[]) 
{
  double h;
  double estimate;
  double sum2;
  double x;
  int i, n;
  int nthreads, tid;
  int procs;
  
  n = atoi(argv[1]);
  h = 1.0 / ( double ) ( 2 * n );
  sum2 = 0.0;

   procs = omp_get_num_procs();
   omp_set_num_threads(procs);  
/* Create a team of threads giving them their own copies of variables */
#pragma omp parallel private( i, x ) shared( h, n ) reduction ( +: sum2 )
{
  printf("h = %lf n = %d\n",h,n);
#pragma omp for
  for ( i = 1; i <= n; i++ )
  {
    x = h * ( double ) ( 2 * i - 1 );
    sum2 = sum2 + 1.0 / ( 1.0 + x * x );
  }
  printf("sum2 is %lf\n", sum2);
}
  estimate = 4.0 * sum2 / ( double ) ( n );
  printf("The estimate of pi is %lf\n",estimate);

  return (0);

}