/*
 * numerical integration example, as discussed in textbook:  
 *
 * compute pi by approximating the area under the curve f(x) = 4 / (1 + x*x)
 * between 0 and 1.
 *
 * parallel version using OpenMP with approach taken in MPI program ("SPMD")
 * and an array of partial sums.  possible performance problems because of
 * "false sharing" of array.
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* copied from not-strictly-standard part of math.h */
#define M_PI		3.14159265358979323846
#include <omp.h>		/* OpenMP header file */

#include "timer.h"

/* main program */
int main(int argc, char *argv[]) {

    char* usage_fmt = "usage:  %s number_of_steps\n";

    /* process command-line arguments */
    char* end_ptr_for_strtol;
    if (argc != 2)  {
        fprintf(stderr, usage_fmt, argv[0]);
        exit(EXIT_FAILURE);
    }
    long num_steps = strtol(argv[1], &end_ptr_for_strtol, 10);
    if (*end_ptr_for_strtol != '\0') {
        fprintf(stderr, usage_fmt, argv[0]);
        exit(EXIT_FAILURE);
    }
 
    double start_time, end_time;
    double pi;
    double sum = 0.0;
    double step = 1.0/(double) num_steps; 
    int nthreads;
    double *partsum;

    /* record start time */
    /* could use omp_get_wtime but get_time is what sequential code uses */
    start_time = get_time();

    /* get number of threads */
    #pragma omp parallel 
    {
        #pragma omp single 
        {
            nthreads = omp_get_num_threads();
        }
    }
    partsum = malloc(sizeof(*partsum) * nthreads);
    if (partsum == NULL) {
        fprintf(stderr, "could not allocate memory\n");
        exit(EXIT_FAILURE);
    }
 
    /* do computation -- using all available threads */
    #pragma omp parallel 
    {
        /* local variables -- or declare in main and make private */
        int myid = omp_get_thread_num();
        double x; 

        partsum[myid] = 0.0;
        for (int i = myid; i < num_steps; i += nthreads) {
            x = (i+0.5)*step;
            partsum[myid] += 4.0/(1.0+x*x);
        }
    }
    /* sum partial results and compute final result */
    for (int k = 0; k < nthreads; ++k) sum += partsum[k];
    pi = step * sum;

    free(partsum);

    /* record end time */
    end_time = get_time();

    /* print results */
    printf("parallel program (SPMD v1a) results with %d threads and %ld steps:\n", 
            nthreads, num_steps);
    printf("computed pi = %g  (%17.15f)\n",pi, pi);
    printf("difference between computed pi and math.h M_PI = %17.15f\n", 
            fabs(pi - M_PI));
    printf("time to compute = %g seconds\n", end_time - start_time);

    return EXIT_SUCCESS;
}