/*
* threadpi.c
* Neil Gershenfeld 12/17/18
* pi calculation benchmark
* pi = 3.14159265358979323846
*/

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

unsigned int npts = 1e8;
double pi = 0;
pthread_mutex_t mutex;

struct data {
   unsigned int index;
   unsigned int points;
   };

void *fn(void *arg) {
   struct data *vars;
   unsigned int i,start,end;
   double sum = 0;
   double a = 0.5;
   double b = 0.75;
   double c = 0.25;
   vars = (struct data *) arg;
   start = 1+(vars->points)*(vars->index);
   end = (vars->points)*(vars->index+1);
   for (i = start; i <= end; ++i)
      sum += a/((i-b)*(i-c));
   pthread_mutex_lock(&mutex);
   pi += sum;
   pthread_mutex_unlock(&mutex);
   pthread_exit(0);
   }

void main(int argc, char *argv[]) {
   int i,nthreads;
   nthreads = atoi(argv[1]);
   pthread_t threads[nthreads];
   struct data var[nthreads];
   double dt,mflops;
   void *status;
   struct timespec tstart,tend;
   clock_gettime(CLOCK_REALTIME,&tstart);
   for (i = 0; i < nthreads; ++i) {
      var[i].index = i;
      var[i].points = npts;
      pthread_create(&threads[i],NULL,fn,(void *) &var[i]);
      }
   for (i = 0; i < nthreads; ++i) {
      pthread_join(threads[i],&status);
      }
   clock_gettime(CLOCK_REALTIME,&tend);
   dt = (tend.tv_sec+tend.tv_nsec/1e9)-(tstart.tv_sec+tstart.tv_nsec/1e9);
   mflops = nthreads*(npts*5.0/(dt*1e6));
   printf("npts = %d, pi = %f\n",npts,pi);
   printf("time = %f, estimated MFlops = %f\n",dt,mflops);
   pthread_mutex_destroy(&mutex);
   pthread_exit(NULL);
   }