Evaluating π using Monte Carlo methods - Serial vs OMP

Question

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;

  #pragma omp parallel private(x, y) num_threads(omp_get_max_threads())
  {

    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?

Answer 1

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 2³¹ - 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.