/*
 * numerical integration example.
 *
 * compute pi by approximating the area under the curve f(x) = 4 / (1 + x*x)
 * between 0 and 1.
 *
 * parallel version using POSIX threads.
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* best(?) available value of pi, to check computed approximation */
/* surprisingly, no library-defined constant for this in standard C! */
#define M_PI acos(-1.0)

#include <pthread.h>

#include "timer.h" 

void * thread_fcn(void * thread_arg);

/* ---- data structure for thread arguments ---- */

typedef struct thread_args {
    int myID;
    int numthreads;
    long numslices;
    double step;
    double *sum_p;
    pthread_mutex_t *lock_p; 
} thread_args_t;

/* ---- main program ---- */

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

    if (argc < 3) {
        fprintf(stderr, "usage %s numthreads\n", argv[0]);
        return EXIT_FAILURE;
    }
    char *endptr;
    long numslices = strtol(argv[1], &endptr, 10);
    if ((*endptr != '\0') || (numslices <= 0)) {
        fprintf(stderr, "numslices must be positive integer\n");
        return EXIT_FAILURE;
    }

    int numthreads = strtol(argv[2], &endptr, 10);
    if ((*endptr != '\0') || (numthreads <= 0)) {
        fprintf(stderr, "numthreads must be positive integer\n");
        return EXIT_FAILURE;
    }

    double start_time, end_time;
    double pi;
    double sum = 0.0;
    double step = 1.0/(double) numslices; 

    thread_args_t thread_args[numthreads];
    pthread_t threads[numthreads];

    pthread_mutex_t lock; 

    /* record start time */
    start_time = get_time();

    /* initialize lock */
    pthread_mutex_init(&lock, NULL);

    /* set up parameters for threads */
    for (int i = 0; i < numthreads; ++i) {
        thread_args_t * args_p = &thread_args[i];
        args_p->myID = i;
        args_p->numthreads = numthreads;
        args_p->numslices = numslices;
        args_p->step = step;
        args_p->sum_p = &sum;
        args_p->lock_p = &lock;
    }

    /* start threads */
    for (int i = 0; i < numthreads; ++i) 
        pthread_create(&threads[i], NULL, thread_fcn, (void *) &thread_args[i]);

    /* wait for all threads to complete */
    for (int i = 0; i < numthreads; ++i) 
        pthread_join(threads[i], NULL);

    /* clean up */
    pthread_mutex_destroy(&lock);

    /* finish computation */
    pi = step * sum;

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

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

    return EXIT_SUCCESS;
}

/* ---- code to be executed by each thread ---- */

void * thread_fcn(void * thread_arg) {
    thread_args_t * args_p = (thread_args_t *) thread_arg;

    double x;
    double part_sum = 0.0;

    /* do this thread's part of computation */
    for (int i=args_p->myID; i < args_p->numslices; i += args_p->numthreads) { 
        x = (i+0.5)*args_p->step;
        part_sum += 4.0/(1.0+x*x); 
    }
    /* add partial sum to global sum */
    pthread_mutex_lock(args_p->lock_p);
    *(args_p->sum_p) += part_sum;
    pthread_mutex_unlock(args_p->lock_p);

    pthread_exit((void* ) NULL);
}