OpenMP parallelisation of pi calculation is either slow or wrong

Question

I'm having trouble parallelising my monte carlo method to calculate pi. Here is the parallelised for-loop:

#pragma omp parallel for private(i,x,y) schedule(static) reduction(+:count)
  for (i = 0; i < points; i++) {
    x = rand()/(RAND_MAX+1.0)*2 - 1.0;
    y = rand()/(RAND_MAX+1.0)*2 - 1.0;

    // Check if point lies in circle
    if(x*x + y*y < 1.0) { count++; }
  }

The problem is, it underestimates pi if I use schedule(static), and its slower than the serial implementation if I use schedule(dynamic). What am I doing wrong? I've tried other ways to fix it (like this: Using OpenMP to calculate the value of PI) but it's still much slower than the serial implementation.

Thanks in advance

Answer 1

Assuming you're using the C library rand function, that function is not reentrant or thread-safe. POSIX provides a rand_r function, but (to quote the glibc documentation):

POSIX.1 extended the C standard functions to support reproducible random numbers in multi-threaded programs. However, the extension is badly designed and unsuitable for serious work.

In particular, the seed must be an unsigned int, which doesn't have enough bits for a good PRNG. They recommend using the SVID random number functions, of which nrand48_r is probably what you're looking for.

Alternatively, you could use a different library.