/* ===================================
   Thread exmaple: compute Pi by integrating the function 1/SQRT(1-x^2)
   =================================== */

#include <iostream.h>
#include <pthread.h>
#include <math.h>

pthread_mutex_t  reduction_mutex;
pthread_t	*tid;

int n, num_threads;
double pi, w;

double f(double a)
{
   return( 2.0 / sqrt(1 - a*a) );
}

void *PIworker(void *arg)
{
   int i, myid;
   double sum, mypi, x;

   /* ----------
      Set individual id to start at 0
      ---------- */
   myid = pthread_self() - tid[0];

   /* ----
      Integrate function f(a), doing only my share !

      Each thread sums the integral usingthe rectangle rule,
      but skips num_threads points each time...
      ---- */
   mypi = 0.0;

   for (i = myid + 1; i <= n; i = i + num_threads)
   {
     x = w * ((double) i - 0.5);
     mypi = mypi + w*f(x);
   }

   /* -----------------
      Thread is done, need to reduce the sum values
      ----------------- */
   pthread_mutex_lock(&reduction_mutex);

   pi = pi + mypi;       /* Add my share to the total... */

   pthread_mutex_unlock(&reduction_mutex);

   return(NULL);     /* Thread exits (dies) */
}

/* =======================
   MAIN
   ======================= */

int main(int argc, char *argv[])
{
   int i;

   /* -----
      Check command line
      ----- */
   if ( argc != 3 )
   {
      cout << "Usage: " << argv[0] << " Num_Intervals Num_Threads" << endl;
      exit(1);
   }

   /* -----
      Get number of intervals and number of threads
      ----- */
   n = atoi(argv[1]);
   num_threads = atoi(argv[2]);


   w = 1.0/(double) n;
   pi = 0.0;

   tid = (pthread_t *) calloc(num_threads, sizeof(pthread_t));


   /* -----
      Initialize the mutex lock
      ----- */
   if ( pthread_mutex_init(&reduction_mutex, NULL) )
   {
      cout << "Cannot init mutex lock" << endl;
      exit(1);
   }

   /* ------
      Create threads
      ------ */
   for (i = 0; i < num_threads; i = i + 1)
   {
      if ( pthread_create(&tid[i], NULL, PIworker, NULL) )
      {
	 cout << "Cannot create thread" << endl;
	 exit(1);
      }
   }

   /* ------
      Once the threads are created, they start compute Pi...
      Each thread will die after they finish its job

      All we need to do is to WAIT for all threads to exit...
      This will make the main thread wait for all thread to exit:
      ------ */
   for (i = 0; i < num_threads; i = i + 1)
      pthread_join(tid[i], NULL);

   /* =======
      The following will only be executed after all therads have exited...
      ======= */

   cout << "Computed Pi = " << pi << endl << endl;

}