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

#define NUM_THREADS 2
#define ITERATIONS 10000000

typedef struct monte_carlo_context_t
{
    int iterations;
    double result;
} mc_context, *pmc_context;

int timespec_subtract (struct timespec *result, struct timespec *x, struct timespec *y)
{
  /* Perform the carry for the later subtraction by updating y. */
  if (x->tv_nsec < y->tv_nsec) {
    int nsec = (y->tv_nsec - x->tv_nsec) / 1000000000 + 1;
    y->tv_nsec -= 1000000000 * nsec;
    y->tv_sec += nsec;
  }
  if (x->tv_nsec - y->tv_nsec > 1000000000) {
    int nsec = (x->tv_nsec - y->tv_nsec) / 1000000000;
    y->tv_nsec += 1000000000 * nsec;
    y->tv_sec -= nsec;
  }

  /* Compute the time remaining to wait.
     tv_usec is certainly positive. */
  result->tv_sec = x->tv_sec - y->tv_sec;
  result->tv_nsec = x->tv_nsec - y->tv_nsec;

  /* Return 1 if result is negative. */
  return x->tv_sec < y->tv_sec;
}

// Generate random numbers and count those inside the circle
void *pi_monte_carlo(void* context)
{
    int i;
    double rx, ry, rs, rt;
    int count = 0;
    pmc_context ctx = (pmc_context)context;
    srand(time(NULL));

    for(i = 0; i < ctx->iterations; i++)
    {
        rx = rand() / (double)RAND_MAX;
        ry = rand() / (double)RAND_MAX;

        double distance = sqrt(rx*rx + ry*ry);
        if (distance <= 1.0)
        {
            count++;
        }
    }

    ctx->result = 4.0 * (count / (double)ctx->iterations);
}

int main(int argc, char* argv[])
{
    int i;
    double result = 0.0;
    struct timespec clk_start, clk_end, clk_diff;

    // Start timer.
    clock_gettime(0, &clk_start);

    // Threads and contexts
    pthread_t threads[NUM_THREADS];
    mc_context contexts[NUM_THREADS];

    for (i = 0; i < NUM_THREADS; i++)
    {
        contexts[i].iterations = ITERATIONS / NUM_THREADS;
        pthread_create(&threads[i], NULL, pi_monte_carlo, (void*)&contexts[i]);
    }

    for (i = 0; i < NUM_THREADS; i++)
    {
        pthread_join(threads[i], NULL);
        result += contexts[i].result;
    }

    // End timer.
    clock_gettime(0, &clk_end);
    timespec_subtract(&clk_diff, &clk_end, &clk_start);

    result /= NUM_THREADS;
    printf("\nThreads: %d, Elapsed time: %d.%09d.\t|\tResult: %lf.\n", NUM_THREADS, clk_diff.tv_sec, clk_diff.tv_nsec, result);
}