# Parallel implementation of „Monte Carlo” algorithm – reactive way

Monte Carlo method lets us to solve some kind of problems using statistic. One of them is to estimate value of PI number (just like trial and error method). To achieve this, we are going to place N random points inside a square (which radius equals 1) and count the number of points that are inside the circle.

### How it works?

There is only one formula for the area of a circle `A = πr2`. Our circle has the radius 1, so it makes that area is just `π`. We can split our circle on 4 parts, so area of each part will be `π/4`. The rest is just a statistic. Take a look at the image below:

This is exactly what we are gonna do. Get N  random points and check how many of them is inside. Then, multiply our result by 4 and voilà! PI number is estimated.

### What we will use?

To solve this problem, we are going to use `windows.h` library, which give us old and really ugly threads API. What is important – our algorithm have to be reactive. Observer pattern should helps us with this. We want to react on „messages” from threads and when all of them will end their calculations – show the result on the screen.

It is not a good solution, critical sections will makes our program slower, we know this. What makes this idea different? The point is that main thread can work in background.

### Let’s start coding

Firstly, we must declare interfaces.

``````
class Observable
{
public:
virtual void notify(const int) = 0;
};

class Observer
{
public:
virtual void onNext(const int) = 0;
virtual void onComplete() const = 0;
};
``````

`Observable` and `Observer` should provide methods, which allows to communicate between them. I think everything is clear here, so now we have to prepare our thread class, which will encapsulate real thread handle.

``````
{
public:

void notify(const int);
void run();

inline const HANDLE getHandle() const { return handle; }

private:
static DWORD WINAPI start(void *);

HANDLE handle;
Observer& observer;
int id count, result;
};
``````

So this is our `ObservableThread`, which contains some properties:

• reference to `Observer` (just one, not list in our case!),
• id – thread universal identifier,
• count – number of iterations,
• result – number of points in circle area.

Okay, now our application class:

``````
class ParallelMCSolver : public Observer
{
public:
ParallelMCSolver(const int, const int);
~ParallelMCSolver();

static const double random();
static const bool isInArea(const double, const double, const double);

void solve();
void waitAll() const;

void onNext(const int);
void onComplete() const;

private:
HANDLE * handlers;
};
``````

`ParallelMCSolver` contains:

• threads vector, which will store our `ObservableThread` instances,
• handlers array (only to make `waitAll()` method work properly),
• nRunning – number of running threads,
• count – number of iterations,
• in – number of points in circle from all threads.

### Uninteresting stuff

We created class definitions, now is time to write a lot of uniteresting code, which do not need explanation (I think), such as constructors, destructors and simple methods.

``````
observer(observer),
handle(NULL),
count(count),
result(0),
id(id)
{}

{
if (handle != NULL) CloseHandle(handle);
}

{
observer.onNext(val);
}

ParallelMCSolver::ParallelMCSolver(const int nThreads, const int count) :
count(count),
in(0)
{
if (nThreads <= 0 || count <= 0) throw exception();

const int tCount = count / nThreads;
for (int i = 0; i < nThreads; i++)
{
}
}

ParallelMCSolver::~ParallelMCSolver()
{
delete[] handlers;
}

void ParallelMCSolver::solve()
{
int i = 0;
{
it->run();
handlers[i] = it->getHandle();
}
}

void ParallelMCSolver::waitAll() const
{
}

const double ParallelMCSolver::random()
{
return (double)rand() / (RAND_MAX + 1);
}

const bool ParallelMCSolver::isInArea(const double x, const double y, const double r)
{
return ((x * x) + (y * y)) <= (r * r);
}

``````

### Make it work

We see that `ParallelMCSolver::solve()` method invokes `run()` of every thread. How `run()` really looks?

``````
{
NULL, 0,
static_cast<void*>(this),
0, NULL);

if (handle == NULL) throw exception();
}
``````

It just create thread using windows API and pass to it `this` pointer of instance, what means, that each thread has access to properties and method of his owner. The start point for thread is `start(void*)` static method.

``````
{
double x, y;

int in = 0;
for (int i = 0; i < thread->count; i++)
{
x = ParallelMCSolver::random();
y = ParallelMCSolver::random();

if (ParallelMCSolver::isInArea(x, y, 1.0)) ++in;
}
return 0;
}
``````

It is the most important function. What’s going on here? Firstly, we have to cast `arg`, which is thread owner – `ObservableThread`.

``````
``````

Then, finally (!) we can take care of our algorithm implementation:

1. Generate random point (x, y).
2. Check if point is inside the circle.
3. If yes – increment `in` variable.
``````
x = ParallelMCSolver::random();
y = ParallelMCSolver::random();

if (ParallelMCSolver::isInArea(x, y, 1.0)) ++in;
``````

After loop, thread has to notify `Observer` and pass him `in` variable. As we know, `notify()` invokes `onNext()` method.

``````
void ParallelMCSolver::onNext(const int value)
{
{
Lock lock(cs);
in += value;
if (--nRunning == 0) onComplete();
}
}
``````

I decided to add some syntatic sugar here. A `Lock` object, which is responsible for enter and leave critical section (implementation below).

``````
class Lock
{
public:
Lock(CRITICAL_SECTION&);
~Lock();

private:
CRITICAL_SECTION& cs;
};

Lock::Lock(CRITICAL_SECTION& cs) : cs(cs)
{
EnterCriticalSection(&cs);
}

Lock::~Lock()
{
LeaveCriticalSection(&cs);
}
``````

When `nRunning` equals 0, which means that all threads finished theirs tasks, `onComplete()` method is called.

``````
void ParallelMCSolver::onComplete() const
{
double PI = 4.0 * (in / (double)count);
cout << "PI: " << PI << endl << endl;
}
``````

`onComplete()` shows us PI result.

### Summary

It’s time to test our code. The last thing is – of course – `main()` function.

``````
CRITICAL_SECTION cs;
int main()
{
const int count = 10000000;

InitializeCriticalSection(&cs);
cout << fixed << setprecision(8);
try
{
while (true)
{

instance.solve();

/*
...
background stuff...
...
*/

instance.waitAll();
}
}
catch (exception& e)
{
cout << e.what() << endl;
}

DeleteCriticalSection(&cs);
cin.get();
return 0;
}
``````

Full source code available at GitHub.