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 &lt; interval; i ++) {
    x = (double)(rand() / (double) RAND_MAX);
    y = (double)(rand() / (double) RAND_MAX);

    if (pow(x,2) + pow(y, 2) &lt;= 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 &lt; interval; i++ ) {
      x = (double)(rand() / (double) RAND_MAX);
      y = (double)(rand() / (double) RAND_MAX);

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

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


This is an efficient method or there is a most efficient version always using OpenMP?
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 &gt; 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 &lt;= 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.