/*
 * 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 POSIX threads.
 *
 * command-line argument gives number of threads.
 */
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>	    /* has pthread_ routines */
#include "timer.h"          /* has get_time() */
#define NUM_STEPS 100000000 

void * thread_fcn(void * thread_arg);

/* ---- global variables (ugly but simple) ---- */
int num_threads;
double sum, step;
pthread_mutex_t lock; 

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

    double start_time, end_time;
    int i;
    int * threadIDs;
    pthread_t * threads;
    double pi;

    num_threads = 0;
    sum = 0.0;
    step = 1.0/(double) NUM_STEPS; 

    if (argc < 2) {
        fprintf(stderr, "usage:  %s numthreads\n", argv[0]);
        return EXIT_FAILURE;
    }
    num_threads = atoi(argv[1]);
    if (num_threads <= 0) {
        fprintf(stderr, "usage:  %s numthreads\n", argv[0]);
        return EXIT_FAILURE;
    }

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

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

    /* set up IDs for threads */
    threadIDs = malloc(num_threads * sizeof(int));
    for (i = 0; i < num_threads; ++i)
        threadIDs[i] = i;

    /* start threads */
    threads = malloc(num_threads * sizeof(pthread_t));
    for (i = 0; i < num_threads; ++i) 
        pthread_create(&threads[i], NULL, thread_fcn, (void *) &threadIDs[i]);

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

    /* clean up */
    pthread_mutex_destroy(&lock);
    free(threadIDs);
    free(threads);

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

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

    /* print results */
    printf("parallel program results with %d threads:\n", num_threads);
    printf("pi = %g  (%17.15f)\n",pi, 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) {
    int myID = * (int *) thread_arg;

    int i;
    double x;
    double part_sum = 0.0;

    /* do this thread's part of computation */
    for (i=myID; i < NUM_STEPS; i += num_threads) { 
        x = (i+0.5)*step;
        part_sum = part_sum + 4.0/(1.0+x*x); 
    }
    pthread_mutex_lock(&lock);  
    sum += part_sum;
    pthread_mutex_unlock(&lock);  

    pthread_exit((void* ) NULL);
}