#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <time.h>

#define PI25DT 3.141592653589793238462643

#define INTERVALS 10000000

int thread_count = 0;
double pi_global = 0;
pthread_mutex_t pi_lock;

void *pi_calc(void *rank);

int main(int argc, char *argv[])
{
    double time2;
    time_t time1 = clock();

    thread_count = strtol(argv[1], NULL, 10);
    pthread_t threads[thread_count];

    for (long i = 0; i < thread_count; i++)
    {
        pthread_create(&(threads[i]), NULL, pi_calc, (void*)i);
    }

    for (long i = 0; i < thread_count; i++)
    {
        pthread_join(threads[i], NULL);
    }

    pthread_mutex_destroy(&pi_lock);

    double dx = 1.0 / (double)INTERVALS;
    double pi = dx * pi_global;

    time2 = (clock() - time1) / (double)CLOCKS_PER_SEC;

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

    return 0;
}

void *pi_calc(void *rank)
{
    long int i, my_rank = (long)rank, intervals = INTERVALS;
    double x, dx, f, sum;
    long remainder = 0;

    if (my_rank == thread_count-1) {
        remainder = intervals %  thread_count;
    }

    long int intervalStart = my_rank*(intervals/thread_count) + (intervals/thread_count);
    long int intervalEnd = my_rank*(intervals/thread_count);


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

    for (i = intervalStart+remainder ; i > intervalEnd; i--)
    {
        x = dx * ((double)(i - 0.5));
        f = 4.0 / (1.0 + x * x);
        sum = sum + f;
    }

    pthread_mutex_lock(&pi_lock);
    pi_global += sum;
    pthread_mutex_unlock(&pi_lock);

    return NULL;
}