// sample pthread program

//      4-Oct-2001      jf nystrom
/*
        A Monte Carlo integration of a
        quarter circle. This will provide
        an estimate for Pi using the rand()
        function.

        Variable input are a "seed" value
        and number of loop iterations.

        23. Jan 2005
        Modified to use user-written random number function
        January 2006
        Modified to use drand and pthreads.
*/

#include <pthread.h>
#include <semaphore.h>
#include <unistd.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <stdlib.h> // needed for srand48_r() and drand48_r()
#define SHARED 1
#define MAX_THREADS     200
#ifndef _REENTRANT
#define _REENTRANT
#endif

 using namespace std;


/*
   Function prototypes.
*/

 void *T_Rand(void *);
 double area_value(double,double);

//      Global variables for use by threads

 pthread_mutex_t vout, nout, add;

 drand48_data status[MAX_THREADS],
                *pts[MAX_THREADS];
 unsigned seed;
 double limit;
 int num_threads = 0;
 double total_over = 0, total_under = 0;

 const double PI = 3.14159265359;

int main(void)
{
cout << MAX_THREADS << "\n";

        // Ask for seed and loop iterations
        cout << "\nenter seed value: "; cin >> seed;
        cout << "\nenter number of loop iterations: "; cin >> limit;

        // Ask for number of threads and check < PTHREAD_THREADS_MAX
        cout << "\nenter number of threads: "; cin >> num_threads;
        if (num_threads > MAX_THREADS)
        {
          cout << "\n ERROR: More than " << MAX_THREADS
                << " threads specified\n.";
          exit(-1);
        }

        // Create thread variables and initialize mutexes
        pthread_t threads[num_threads];

        pthread_mutex_init(&vout,NULL);
        pthread_mutex_init(&nout,NULL);
        pthread_mutex_init(&add,NULL);

        // Create threads
        int check;

        cout <<  "Creating " << num_threads << " threads...\n";

        for (int i=0; i < num_threads; i++)
        {
          check = pthread_create(&threads[i], NULL, T_Rand, (void *)(i+1));
          if (check)
          {
            cout << "ERROR; return code from pthread_create() is "
                  << check << "\n";
            exit(-1);
          }
          else
          {
            pthread_mutex_lock(&nout);
            cout << "Thread " << i+1 << " created.\n";
            pthread_mutex_unlock(&nout);
          }
        }

        // Join threads after termination
        for (int i=0; i < num_threads; i++)
          pthread_join(threads[i], NULL);

        // Print totals
        double my_value = 4.*area_value(total_over,total_under),
           my_error = sqrt((PI-my_value)*(PI-my_value));

        cout << "\n\nThe value (using " << num_threads
          << " threads) of Pi = " << setprecision(10) << fixed
          << my_value << "\n  with error of " << my_error
          << "\n  (" << setprecision(0) << fixed
          << limit << " iterations)" << "\n\n";

}
void *T_Rand(void *t)
{
        int worker = (int) t;
        double over = 0, under = 0;
        double x, y, *ptx, *pty;;

        pthread_mutex_lock(&nout);
        cout << "\tThread " << worker << " started.\n";
        pthread_mutex_unlock(&nout);

        // seed local generator
        pts[worker] = &status[worker];
        srand48_r(seed*worker,pts[worker]);

        // count number of points above and below for this thread

        ptx = &x; pty = &y;
        double count = 0;

        for (double n=worker; n <= limit; n+= num_threads)
        {
          count++;
          drand48_r(pts[worker],ptx);
          drand48_r(pts[worker],pty);

          if (y > sqrt(1. - x*x)) over++;
          else under++;
        }

        // Print local values. (Need to mutex for use of function.)

        pthread_mutex_lock(&vout);

        double my_value = 4.*area_value(over,under),
           my_error = sqrt((PI-my_value)*(PI-my_value));

        pthread_mutex_unlock(&vout);

        pthread_mutex_lock(&nout);

        cout << "\n\nThe value (for thread " << worker
          << " ) of Pi = " << setprecision(10) << fixed
          << my_value << "\n  with error of " << my_error
          << "\n  ( " << setprecision(0) << fixed
          << count << " iterations)" << "\n\n";

        pthread_mutex_unlock(&nout);

        // Add local totals to global totals

        pthread_mutex_lock(&add);

        total_over += over;
        total_under += under;

        pthread_mutex_unlock(&add);

        pthread_exit(NULL);

}


/*
        This function calculates the Monte Carlo area
*/

double area_value(double over, double under)
{
        double value;
        value = under/(over + under);
        return value;
}

/*

$g++ -O -pthread monte.cpp
$a.out

enter seed value: 345

enter number of loop iterations: 1000000000

enter number of threads: 4
Creating 4 threads...
Thread 1 created.
        Thread 1 started.
Thread 2 created.
        Thread 2 started.
Thread 3 created.
        Thread 3 started.
Thread 4 created.
        Thread 4 started.

The value (for thread 1 ) of Pi = 3.141508448000
  with error of 0.000084205590
  (250000000 iterations)



The value (for thread 4 ) of Pi = 3.141626928000
  with error of 0.000034274410
  (250000000 iterations)



The value (for thread 3 ) of Pi = 3.141569120000
  with error of 0.000023533590
  (250000000 iterations)



The value (for thread 2 ) of Pi = 3.141716544000
  with error of 0.000123890410
  (250000000 iterations)



The value (using all threads) of Pi = 3.141605260000
  with error of 0.000012606410
  (1000000000 iterations)


*/