Search code examples
crandomopenmprandom-seed

How to generate uniformly distributed random numbers between 0 and 1 in a C code using OpenMP?


I am trying to write an OpenMP code in which each thread will work on big arrays of uniformly distributed random numbers between 0 and 1. Each thread needs to have different and independent random number distributions. In addition, the random number distributions need to be different every time the code is called. This is what I am using right now. Does this always guarantee each thread has its own/different random number sequences? Will the sequences be different every time the code is called? What is the correct way of doing this? The following code has each thread generating 5 samples but in an actual run it will be order of millions.

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <time.h>

int main(int argc, char* argv[])
{
    int numthreads,i;
    #pragma omp parallel private(i)
    {
        int id;
        id=omp_get_thread_num();
        if(id==0) numthreads = omp_get_num_threads();
        printf("thread %d \n",id);
        srand(time(0)^omp_get_thread_num());
        for (i=0; i<5; i++)
        {
            printf("thread %d: %d %.6f \n",id,i,(double)rand()/(double)RAND_MAX);
        }
    }
    return 0;
}

Solution

  • You don't mention what OS you're using, but if it's Linux or a POSIX compliant system, there's erand48() for thread-safe generation of random numbers uniformly distributed in the range [0.0, 1.0). It uses a 48-bit seed that's passed as an argument. Generating the initial seed can be done in a number of ways. OpenBSD and Linux have getentropy(), BSDs have arc4random_buf(), you can read from the /dev/urandom special file on many OSes, or do something like you're currently using with time, pid, etc. I'd suggest a higher resolution timer than time(), though - clock_gettime() is a good source.

    An example:

    #include <errno.h>
    #include <fcntl.h>
    #include <omp.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include <sys/stat.h>
    #include <sys/types.h>
    #include <unistd.h>
    
    int main(void) {
    #pragma omp parallel for
      for (int i = 0; i < 4; i++) {
        unsigned short xi[3]; // PRNG state variable
    
    #if 0
        // OpenBSD 5.6+/Linux kernel 3.17+ and glibc 2.25+
        if (getentropy(xi, sizeof xi) < 0) {
          perror("getentropy");
          exit(EXIT_FAILURE);
        }
    #else
        // Read from /dev/urandom
        int fd = open("/dev/urandom", O_RDONLY);
        if (fd < 0) {
          perror("open /dev/urandom");
          exit(EXIT_FAILURE);
        }
        if (read(fd, xi, sizeof xi) != sizeof xi) {
          perror("read");
          exit(EXIT_FAILURE);
        }
        close(fd);
    #endif
    
        for (int n = 0; n < 4; n++) {
          printf("Thread %d random number %f\n", omp_get_thread_num(), erand48(xi));
        }
      }
    
      return 0;
    }