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

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

#include <sys/time.h>

struct timeval start_time, stop_time;
int    elapsed;

pthread_mutex_t  sum_mutex;

int N, num_threads;
double sum, w;

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

void *PIworker(void *arg)
{
   int i, myStart;
   double x, tmp;

   /* ----------
      Get starting index
      ---------- */
   myStart = * (int *) arg;

   /* ----
      Add all partial sum to tmp before adding to sum !!!
      ---- */
   tmp = 0.0;
   for (i = myStart; i < N; i = i + num_threads)
   {
     x = w * ((double) i + 0.5);
     tmp = tmp + w*f(x);
   }

   pthread_mutex_lock(&sum_mutex);
   sum = sum + tmp;			// Synchr only ONCE per thread
   pthread_mutex_unlock(&sum_mutex);

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

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

int main(int argc, char *argv[])
{
   int i;
   pthread_t tid[100];
   int Start[100];       // Start index values

   /* -----
      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;
   sum = 0.0;

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


   gettimeofday(&start_time, NULL);		// Get time

   /* ------
      Create threads
      ------ */
   for (i = 0; i < num_threads; i = i + 1)
   {
      Start[i] = i;
      if ( pthread_create(&tid[i], NULL, PIworker, &Start[i]) )
      {
	 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);


   gettimeofday(&stop_time, NULL);		// Get time

   elapsed = (stop_time.tv_sec*1000000 + stop_time.tv_usec) -
                (start_time.tv_sec*1000000 + start_time.tv_usec);


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

   cout << "Computed Pi = " << sum << endl << endl;
   cout << "Elapsed time = " << elapsed << " microseconds" << endl;

}