Solution 3

C

/* Exercise: Pi
 *
 * In this exercise you will determine the value 
 * of PI using the integral  of 
 *    4/(1+x*x) between 0 and 1. 
 *
 * The integral is approximated by a sum of n intervals.
 * 
 * The approximation to the integral in each interval is:
 *    (1/n)*4/(1+x*x). 
 */
#include <stdio.h>
#include <time.h>
#ifdef _OPENMP
#include <omp.h>
#endif

#define PI25DT 3.141592653589793238462643
#define INTERVALS 100000000

int main(int argc, char **argv)
{
    long int i, intervals = INTERVALS;
    double x, dx, f, sum, pi;
    double time2;
#ifdef _OPENMP 
    double time1 = omp_get_wtime();
#else
    time_t time1 = clock();
#endif
    
    printf("Number of intervals: %ld\n", intervals);

    sum = 0.0;
    dx = 1.0 / (double) intervals;

    #pragma omp parallel for private(x,f) \
    reduction(+:sum)
    for (i = 1; i <= intervals; i++) {
        x = dx * ((double) (i - 0.5));
        f = 4.0 / (1.0 + x*x);
        sum = sum + f;
    }
         
    pi = dx*sum;

#ifdef _OPENMP
    time2 = omp_get_wtime() - time1;
#else
    time2 = (clock() - time1) / (double) CLOCKS_PER_SEC;
#endif

    printf("Computed PI %.24f\n", pi);
    printf("The true PI %.24f\n\n", PI25DT);
    printf("Elapsed time (s) = %.2lf\n", time2);

    return 0;
}

FORTRAN

!---------------------------------------------c
!  Exercise: Pi                               c
!                                             c
!  Compute the value of PI using the integral c
!  pi = 4* int 1/(1+x*x)    x in [0-1]        c
!                                             c
!  The integral is approximated by a sum of   c
!  n interval.                                c
!                                             c
!  The approximation to the integral in each  c
!  interval is: (1/n)*4/(1+x*x).              c
!---------------------------------------------c
program pigreco
#ifdef _OPENMP
    use omp_lib
#endif
    implicit none

    integer(selected_int_kind(18)) :: i 
    integer(selected_int_kind(18)), parameter :: intervals=1e8

    real(kind(1.d0)) :: dx,sum,x
    real(kind(1.d0)) :: f,pi

    real(kind(1.d0)), parameter :: PI25DT = acos(-1.d0)

#ifdef _OPENMP
    real(kind(1.d0)) :: time1, time2
    time1 = omp_get_wtime()
#else 
    real :: time1, time2
    call cpu_time(time1)
#endif

    print *, 'Number of intervals: ', intervals
    sum=0.d0
    dx=1.d0/intervals

!$omp parallel do private(x,f) &
!$omp reduction(+:sum)
    do i=1,intervals
        x=dx*(i-0.5d0)
        f=4.d0/(1.d0+x*x)
        sum=sum+f
    end do
!$omp end parallel do

    pi=dx*sum
#ifdef _OPENMP
    time2 = omp_get_wtime()
#else
    call cpu_time(time2)
#endif

    PRINT '(a13,2x,f30.25)',' Computed PI =', pi
    PRINT '(a13,2x,f30.25)',' The True PI =', PI25DT
    PRINT *, ' '
    PRINT *, 'Elapsed time ', time2-time1 ,' s' 

end program