4
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ Are you asking whether the Monte Carlo approach is efficient or the implementation with OpenMP? \$\endgroup\$ – Ninetails Jul 14 at 21:01
  • \$\begingroup\$ @Ninetails My implementation in OpenMp \$\endgroup\$ – Ilmionome456 Jul 14 at 23:16
1
\$\begingroup\$

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.
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.