# Evaluating π using Monte Carlo methods - Serial vs OMP

I wrote this simple code for evaluating the π using Monte Carlo method. This is the serial version:

long double compute_pi_serial(const long interval) {

srand(time(NULL));

double x, y;
int i, circle = 0;

for (i = 0; i < interval; i ++) {
x = (double)(rand() / (double) RAND_MAX);
y = (double)(rand() / (double) RAND_MAX);

if (pow(x,2) + pow(y, 2) <= 1.0) circle ++;
}

return (long double) circle / interval * 4.0;
}


After, I wrote a parallel version using OpenMP and this is the result:

long double compute_pi_omp(const long interval) {

double x, y;
int i, circle = 0;

{

srand(SEED);

#pragma omp for reduction(+:circle)
for ( i = 0; i < interval; i++ ) {
x = (double)(rand() / (double) RAND_MAX);
y = (double)(rand() / (double) RAND_MAX);

if ( pow(x, 2) + pow(y, 2) <= 1 ) circle++;
}
}

return (long double) circle / interval * 4.0;
}


This is an efficient method or there is a most efficient version always using OpenMP?

• Are you asking whether the Monte Carlo approach is efficient or the implementation with OpenMP? Jul 14, 2019 at 21:01
• @Ninetails My implementation in OpenMp Jul 14, 2019 at 23:16

Correctness issues

Correctness merits priority over efficiency. Important to get a correct answer, not a fast and wrong one.

• Poor use of mixing types

long interval, and int i lead to undefined behavior (UB) when interval > INT_MAX. Use same type.

• Small range

long may use 32-bit. Easy enough to run in less than 1 minute and exceed 231 - limiting how far we can test this code. Use long long or, use unsigned long long.

• Off-by-1

Using / (double) RAND_MAX may impart a small bias. Yet given the slow convergence of the Monte Carlo PI, I doubt any serious miscalculation.

I'd expect the random values to be in the middle of each random range:

[0.5, 1.5, 2.5, ... RAND_MAX - 1.5, RAND_MAX - 0.5, RAND_MAX + 0.5]/(RAND_MAX + 1u)


Code has then distributed differently:

[0.0, 1.0/RAND_MAX, 2.0/RAND_MAX, ... 1.0]


A more significant bias occurs when (int)(double) RAND_MAX != RAND_MAX - select platforms with more precision in RAND_MAX than double.

Note: to convert RAND_MAX plus 1 to double, code can use (RAND_MAX/2 + 1)*2.0. Works well when RAND_MAX is a Mersenne Number or and odd value. This avoids losing precision that (double)RAND_MAX + 1.0 may incur.

Efficiency

• pow() vs x*x

A weak compiler may use a laborious function call seeing pow(x,2) and not the certainly more efficient x*x replacement. Of course the right answer is to use a smart compiler. Still, profiling a direct coding of x*x may prove worth-wild.

• fp vs. integer

I'd even consider an integer only approach. Something like the below.

// int2x twice as wide as int
int2x limit = (int2x) RAND_MAX * RAND_MAX;
...
int2x x = rand();
int2x y = rand();
if (x*x + y*y <= limit) circle++;

• double vs. long double

Seriously doubt any calculation run within a day will benefit using the higher precision long double.

// return (long double) circle / interval * 4.0;
return circle / interval * 4.0;

• No comment on OP's major goal: Serial vs OMP - sorry.