/* File:       pi-pthread.c
 *
 * Purpose:    Estimates pi in parallel using the Madhava–Leibniz series and
 *             pthreads. NOTE: the number of intervals should be
 *             evenly-divisible by the number of threads.
 *
 * Compile:    gcc -g -Wall -pthread -o pi-pthread pi-pthread.c
 *                 (add -DDEBUG to see thread info)
 *                 (bonus: try with -O3 ?)
 * Run:        ./pi-pthread intervals [num-threads]
 */

#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <unistd.h>

/* We define pi here so we can check and see how accurate our computation is. */
#define PI 3.141592653589793238462643

/** Size of the calculation performed by each thread: */
unsigned long long block_size;

long num_threads = 1;
long double *totals;

double get_time(void)
{
    struct timeval t;
    gettimeofday(&t, NULL);
    return t.tv_sec + t.tv_usec / 1000000.0;
}

void *worker_thread(void *arg)
{
    unsigned long rank = (unsigned long) arg;
    unsigned long long start = block_size * rank;
    unsigned long long end = block_size * rank + block_size;
#ifdef DEBUG
    fprintf(stderr, "I am thread %ld calculating range %Ld to %Ld\n",
            (long) rank, start, end);
#endif

    for (unsigned long long i = start; i < end; ++i) {
        if (i % 2 == 0) {
            totals[rank] = totals[rank] + (1.0 / ((2.0 * i) + 1));
        } else {
            totals[rank] = totals[rank] - (1.0 / ((2.0 * i) + 1));
        }
    }

    return NULL;
}

int main(int argc, char *argv[])
{
    if (argc <= 1) {
        fprintf(stderr, "Usage: %s num-intervals [num-threads]\n", argv[0]);
        return EXIT_FAILURE;
    }
    
    if (argc < 3) {
        num_threads = sysconf(_SC_NPROCESSORS_ONLN);
    } else if (argc >= 3) {
        num_threads = strtol(argv[2], NULL, 10);
    }

    unsigned long long intervals = strtoull(argv[1], NULL, 10);
    printf("Intervals: %llu\n", intervals);
    printf("  Threads: %ld\n", num_threads);

    if (intervals == 0 || num_threads <= 0) {
        return 0;
    }

    block_size = intervals / num_threads;

    double time1 = get_time();

    totals = malloc(num_threads * sizeof(long double));
    pthread_t *threads = malloc(num_threads * sizeof(pthread_t));
    for (int i = 0; i < num_threads; ++i) {
        pthread_create(
                &threads[i],
                NULL,
                worker_thread,
                ((void *) (unsigned long) i));
    }

    long double total = 0.0;
    for (int i = 0; i < num_threads; ++i) {
        pthread_join(threads[i], NULL);
        total += totals[i];
    }
    total = total * 4;

    double time2 = get_time();
    printf("   Result: %.20Lf\n", total);
    printf(" Accuracy: %.20Lf\n", PI - total);
    printf("     Time: %.10lf\n", time2 - time1);
}