Monte Carlo Integration method to calculate PI value (C++/MT)

This is a simple program to calculate the PI value using the new C++11 random facility and MultiThread. Change the number of wanted thread to fit your CPU power.

The code is not pretty well written (threads work forever, output a little messed up), but it runs and could be a good starting point for experimenting.

To know more about Monte Carlo Integration statistical, non deterministic method, read this Wiki article.

/**
Use the MonteCarlo method to calculate the pi-greco value.

*/

#include <iostream>
#include <iomanip>		// setprecision
#define _USE_MATH_DEFINES
#include <math.h>
#include <thread>
#include <mutex>


#if defined(_MSC_VER)
#	include <windows.h>
#	include <time.h>
#else
#	include <sys/timeb.h>
#endif

// New C++11 random generator API
#include <random>

#if !defined(M_PI)
#	define M_PI          3.141592653589793238462643383279502884L /* pi */
#endif

using namespace std;

// Spiegazione del metodo monte-carlo
// http://www.eveandersson.com/pi/monte-carlo-circle
// Libro sul metodo monte carlo e le sue diverse applicazioni nella simulazione
// http://www.javaquant.net/books/MCBook-1.2.pdf
// Utilizzo della nuova funzionalità random in C++11:
// http://stackoverflow.com/questions/19665818/best-way-to-generate-random-numbers-using-c11-random-library

class IRandomDevice
{
public:
	virtual void randomize() = 0;
	virtual double get_data() = 0;
	double operator() ()
	{
		return get_data();
	}
};

typedef long long tCounter;

/// Old C rand/randomize functions
class random_device_old : public IRandomDevice
{
public:
	void randomize()
	{
#if defined(_MSC_VER)
		FILETIME fileStart;
		GetSystemTimeAsFileTime(&fileStart);
		ULARGE_INTEGER start;
		start.HighPart = fileStart.dwHighDateTime;
		start.LowPart = fileStart.dwLowDateTime;
		srand(start.LowPart*start.HighPart);
#else
		struct timeval t1;
		gettimeofday(&t1, NULL);
		srand(t1.tv_usec * t1.tv_sec);
#endif

	}
	double get_data()
	{
		return (double)rand() / RAND_MAX;
	};
};

/// New std C++11 random facility
template <typename RAND_DEVICE = std::random_device, typename RAND_GENERATOR = std::mt19937, typename RAND_DISTRIBUTION = std::uniform_real_distribution<double>>
class random_device_c11 : public IRandomDevice
{
public:
	RAND_DEVICE			rd;
	RAND_GENERATOR		mt;
	RAND_DISTRIBUTION	dist;
	random_device_c11() : mt(rd()), dist(0, 1)
	{

	}
	void randomize() 
	{

	}
	double get_data()
	{
		
		return dist(mt);
	}
};

class CAcc
{
public:
	virtual void accumulate(const char * pTitle, tCounter tot, tCounter in_circle)
	{
		double new_pi = ((double)in_circle / tot)*4.0;
		std::cout << pTitle << ": " << tot/1000000 << "M/" << in_circle << " value= " << new_pi << " D= " << abs(new_pi - M_PI) << std::endl;
	}
};

class CMTAcc : public CAcc
{
public:
	std::mutex m_Mutex;
	tCounter m_tot;
	tCounter m_in;
	CMTAcc()
	{
		m_tot = m_in = 0;
	}
	virtual void accumulate(const char * pTitle, tCounter tot, tCounter in_circle)
	{
		tCounter temp_tot, temp_in_circ = 0;
		{
			std::lock_guard<std::mutex> lk(m_Mutex);
			temp_tot = m_tot += tot;
			temp_in_circ = m_in += in_circle;
		}
		CAcc::accumulate(pTitle, temp_tot, temp_in_circ);
	}
};



void _PI_MonteCarlo(const char *  title, IRandomDevice * rand_gen, CAcc & Acc)
{
	printf("Reference value : %3.15lf\n", M_PI);

	(*rand_gen).randomize();
	tCounter nCounts = 0;
	tCounter n = 0;
	tCounter dn = 10000000;
	
	while (1)
	{
		++n;
		double x = (*rand_gen)();
		double y = (*rand_gen)();
		if (sqrt(x*x + y*y) < (double)1.0) 
			++nCounts;
		if (n%dn == 0)
		{
			Acc.accumulate(title, n, nCounts);
			(*rand_gen).randomize();
		}
	}
	delete rand_gen;
}



std::vector<thread*> joins;
const char * title1 = "rand_method";
const char * title2 = "c++11_method";
void PI_MonteCarlo()
{
	cout << setprecision(14);

	CAcc Acc1;
	int n_threads = 8;
	auto f1 = std::bind(_PI_MonteCarlo, title1, new random_device_old(), std::ref(Acc1));
	std::thread * th1= new std::thread (f1);
	joins.push_back(th1);


	CMTAcc Acc2;
	for (int i = 0; i < n_threads; i++)
	{
		auto f2 = std::bind(_PI_MonteCarlo, title2, new random_device_c11<>(), std::ref(Acc2));
		std::thread * th2= new thread(f2);
		joins.push_back(th2);
	}

	for (auto i = begin(joins); i != end(joins);i++)
	{
		(*i)->join();
	}

}