/*
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
#include
#include
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);
}