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

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

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

      Each thread sums the integral usingthe rectangle rule,
      but skips num_threads points each time...
      ---- */
   for (i = myStart; i < N; i = i + num_threads)
   {
     x = w * ((double) i + 0.5);

     pthread_mutex_lock(&sum_mutex);
     sum = sum + w*f(x);
     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 start 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;

}