# Calculating pi with Monte Carlo using OpenMP

I'm a C noob and I'm learning about concurrency using C. I came across an exercise in a book asking me to find the approximate value of Pi using the Monte Carlo technique with OpenMP. I came up with the following:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

int in_circle = 0;
int total = 10000;

void calculate_in_circle_random_count() {
srand((unsigned int)time(NULL));
for (int i = 0; i < total / thread_numb; i++) {
float x = (float)rand()/(float)(RAND_MAX);
float y = (float)rand()/(float)(RAND_MAX);

float val = (x * x) + (y * y);

if (val < 1.0) {
#pragma omp critical
{
in_circle++;
}
}
}
}

int main(int argc, char const *argv[]) {
{
calculate_in_circle_random_count();
}

float pi_approx = 4 * (float)in_circle / (float)total;

printf("In circle: %d\n", in_circle);
printf("Total: %d\n", total);
printf("Pi approximation: %.6f\n", pi_approx);

return 0;
}


I have the sample count of 10000 to calculate approximation and I want to dedicate 4 threads to be parallel for the calculation. So since I want 4 threads, I make sure the for loop is run until total / 4 for each thread. And whenever (val < 1.0) is true I increase the in_circle variable in a critical section.

Does this approach make sense? If not, how could it be improved?

This code has a serious problem due to (at least the typical implementation of) rand(). rand() normally has a (hidden) seed value, and each call to rand() modifies the seed. At least in most implementations, that means rand always forces serialization. That means, calling rand() (a couple of times) in your inner loop will prevent code from scaling well at all. In fact, in most cases (as Toby Speight showed) a multi-threaded version will run substantially slower than a single-threaded version.

To fix this, you pretty much need to use some other random number generator. Your implementation may provide one (e.g., erand48 is fairly common). If you really need your code to be portable, you could write your own, something on this order:

#include <time.h>

typedef unsigned long long rand_state;

// multiplier/modulus taken from Knuth Volume 2, Table 1
static const int multiplier = 314159269;
static const int addend = 1;
static const int modulus = 0xffffffff;

// note that this works differently from srand, returning a seed rather than setting
// a hidden seed.
rand_state omp_srand() {
rand_state state = time(NULL);
state ^= (unsigned long long)omp_get_thread_num() << 32;
return state;
}

int omp_rand(rand_state *state) {
*state = *state * multiplier + addend;
return *state & modulus;
}


Note that since it's entirely likely that all the threads get started in the same second, this combines the current time with the thread ID to seed each thread's generator. Although it's not 100% guaranteed, this gives a pretty high likelihood that each thread's generator will start with a unique seed. On the other hand, it also means that we have to start the thread, and then seed the generator inside that thread.

To use this, our code would be something on this general order:

double calculate_in_circle_random_count(void) {

static const unsigned total = 1000000000;
unsigned in_circle = 0;

#pragma omp parallel reduction(+:in_circle)
{
rand_state seed = omp_srand();
int count = total / omp_get_num_threads();
int i;

// then do this thread's portion of the computation:
for (i = 0;  i < count;  ++i) {

double x = (double)omp_rand(&seed) / OMP_RAND_MAX;
double y = (double)omp_rand(&seed) / OMP_RAND_MAX;

double val = x * x + y * y;
in_circle +=  val < 1.0;
}
}
return 4.0 * in_circle / total;
}

int main(int argc, char const *argv[])
{

if (argc > 1)

float pi_approx = calculate_in_circle_random_count();

printf("Pi approximation: %.6f\n", pi_approx);

return 0;
}


With this modification, the code is at least capable of scaling. To get it to scale very well, you'd need to change your total to a rather larger number though--with it as small as you've specified, it takes longer to start up multiple threads than it saves in calculation. But at least this code can scale well, so if we make total quite a bit larger, it really will run faster. For timing, I added a few more zeros to total, and rewrote main a bit, to create a loop to use 1, 2, 3, and 4 threads and print out the time each iteration. On my machine, this produced the following:

With 1 threads: Pi approximation: 3.140259, time: 3,652,683 microseconds
With 2 threads: Pi approximation: 3.139389, time: 1,892,415 microseconds
With 3 threads: Pi approximation: 3.138885, time: 1,268,917 microseconds
With 4 threads: Pi approximation: 3.138306, time:   935,579 microseconds


So with this, 4 threads is at least close to 4 times as fast as one thread. I'm pretty sure if you use rand inside the loop, you'll never get it to scale well at all.

For anybody who cares, the version of main that produces that output is written in C++, and looks like this:

int main(int argc, char const* argv[]) {
using namespace std::chrono;

std::cout.imbue(std::locale(""));

int processors = omp_get_num_procs();

std::cout << "With " << std::setw(2) << num_threads << " threads: ";

auto start = high_resolution_clock::now();
float pi_approx = calculate_in_circle_random_count();
auto stop = high_resolution_clock::now();

printf(" %.6f, time: ", pi_approx);
std::cout << std::setw(9) << duration_cast<microseconds>(stop - start).count() << " microseconds\n";
}
return 0;
}


This remains pretty much the same, but with some timing code hacked in. In particular, it's not an attempt at rewriting the code in C++ in general. If I were doing that, I'd probably do a number of things rather differently (starting with the fact that the C++ <random> header already provides a clean way to handle per-thread random number generation, so I'd use that instead of rolling my own).

• Won't this end up with all threads using the same random number seed, because seed is initialized via omp_srand() before the OMP loop is entered, so that initial value will be copied to all the threads? – 1201ProgramAlarm Feb 21 at 3:47
• This is very informative! Would you mind sharing the benchmarking code where you got the time comparisons? I would like to compare how using rand() affects performance locally. – emrepun Feb 21 at 7:07
• @emrepun: I don't mind sharing it, but (at least right now) it's written in C++, because it uses the C++ chrono library for timing. – Jerry Coffin Feb 21 at 7:09
• @emrepun: For what it's worth, I believe OpenMP provides a timing function as well. Oh, and I've edited the version of main with the chrono-based timing code into the answer above. – Jerry Coffin Feb 21 at 7:26
• @emrepun: oops--a copy and paste error on my part. It should be: #define OMP_RAND_MAX 0xffffffff. Just to be clear: it's something I made up for the omp_rand/omp_srand I wrote for this, not something defined by OpenMP. – Jerry Coffin Feb 21 at 7:32

I think you're over-complicating this. You can let OpenMP divide the work between the threads for you, using #pragma openmp parallel for.

You can use a reduction for in_circle instead of a critical section - that allows OpenMP to sum into a per-thread variable, and add them all at the end, reducing contention for the variable.

Here's a simpler version (I've removed the unused includes, too, and made the function declaration a prototype):

#include <stdio.h>
#include <stdlib.h>

double calculate_in_circle_random_count(void)
{
static const unsigned total = 1000000;
unsigned in_circle = 0;

#pragma omp parallel for reduction(+:in_circle)
for (unsigned i = 0;  i < total;  ++i) {
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;

double val = x * x + y * y;
in_circle +=  val < 1.0;
}

return 4.0 * in_circle / total;
}

int main(void)
{
const double pi_approx = calculate_in_circle_random_count();
printf("Pi approximation: %.6f\n", pi_approx);
}


Note that at this scale, the parallelisation overheads dominate, making the parallel version much slower than the serial. With more iterations, you'll get closer, but there's so little in the loop that the benefits are small.

• I see this version is much simpler! One question though, I see you removed srand((unsigned int)time(NULL)); I thought that was needed to properly seed the random number generation. Is it not only needed but would also harm the performance of the code? – emrepun Feb 21 at 7:04
• Do you really need a different result every time? If so, then it's appropriate to call srand(), but here it seems that starting with the generator in its default state is adequate. BTW, I didn't check whether rand() generates distinct values in each thread, so you may find that all the threads are working with the same inputs. – Toby Speight Feb 21 at 8:34
• Worse that that, I have actually checked and discovered rand() is not even thread-safe, so calling it concurrently could actually corrupt its state. You'll need to find a thread-safe generator, such as POSIX rand_r() or its successor erand48() (which conveniently yields values between 0.0 and 1.0). – Toby Speight Feb 21 at 8:36
• There are random number generators which are explicitly designed for use in parallel environments (and are available in many math libraries). Read "Parallel Random Numbers: As Easy as 1, 2, 3 - The Salmons", fhttp://www.thesalmons.org/john/random123/papers/random123sc11.pdf for instance. Merely initialising each thread's generator with a random seed is could still lead to sampling the same sequence in multiple threads if you are unlucky... – Jim Cownie Feb 22 at 8:41