How to Generate a Random Double Uniformly Distributed Between 0 and 1 from C++

How do you generate a random double uniformly distributed between 0 and 1 from C++?

In C++11 and C++14 we have much better options with the random header. The presentation rand() Considered Harmful by Stephan T. Lavavej explains why we should eschew the use of rand() in C++ in favor of the random header and N3924: Discouraging rand() in C++14 further reinforces this point.

The example below is a modified version of the sample code on the cppreference site and uses the std::mersenne_twister_engine engine and the std::uniform_real_distribution which generates numbers in the [0,1) range (see it live):

#include <iostream>
#include <iomanip>
#include <map>
#include <random>

int main()
{
    std::random_device rd;

    std::mt19937 e2(rd());

    std::uniform_real_distribution<> dist(0, 1);

    std::map<int, int> hist;
    for (int n = 0; n < 10000; ++n) {
        ++hist[std::round(dist(e2))];
    }

    for (auto p : hist) {
        std::cout << std::fixed << std::setprecision(1) << std::setw(2)
                  << p.first << ' ' << std::string(p.second/200, '*') << '\n';
    }
}

output will be similar to the following:

0 ************************
1 *************************

Since the post mentioned that speed was important then we should consider the cppreference section that describes the different random number engines (emphasis mine):

The choice of which engine to use involves a number of tradeoffs*: the
**linear congruential engine is moderately fast and has a very small
storage requirement for state. The lagged Fibonacci generators are
very fast even on processors without advanced arithmetic instruction
sets, at the expense of greater state storage and sometimes less
desirable spectral characteristics. The Mersenne twister is slower and
has greater state storage requirements but with the right parameters
has the longest non-repeating sequence with the most desirable
spectral characteristics (for a given definition of desirable).

So if there is a desire for a faster generator perhaps ranlux24_base or ranlux48_base are better choices over mt19937.

rand()

If you forced to use rand() then the C FAQ for a guide on How can I generate floating-point random numbers?, gives us an example similar to this for generating an on the interval [0,1):

#include <stdlib.h>

double randZeroToOne()
{
    return rand() / (RAND_MAX + 1.);
}

and to generate a random number in the range from [M,N):

double randMToN(double M, double N)
{
    return M + (rand() / ( RAND_MAX / (N-M) ) ) ;  
}

Generating random double between 1 and 100 in C

What you're doing would "work", but probably not in the way you intend to. The key is to define what you mean by random double. Normally, what we'd mean by this is a double with the maximal precision the system allows.

Fortunately, we know that the number returned by rand is between 0 and RAND_MAX, whatever that value happens to be. (32767 used to be a common choice.) We can divide by this limit to generate a value between 0 and 1, which then can be scaled to whatever range you need, like so:

double random = ((double) rand() / RAND_MAX) * (100 - 1) + 1;

If you wanted a value between 1.5 and 9, you'd do it in a very similar way:

double random = ((double) rand() / RAND_MAX) * (9 - 1.5) + 1.5;

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

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;
}

Generating a random, uniformly distributed real number in C

You can use:

#include <time.h>
srand(time(NULL)); // randomize seed
double get_random() { return (double)rand() / (double)RAND_MAX; }
n = get_random();

srand() sets the seed which is used by rand to generate pseudo-random numbers. If you don't call srand before your first call to rand, it's as if you had called srand(1) (serves as a default).

If you want to exclude [1] use:

(double)rand() / (double)((unsigned)RAND_MAX + 1);

Full solution:

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

double get_random() { return ((double)rand() / (double)RAND_MAX); }

int main()
{
    double n = 0;
    srand(time(NULL)); // randomize seed
    n = get_random(); // call the function to get a different value of n every time

    printf("%f\n", n);  // print your number
    return 0;
}

Every time you run it you will get a different number for n.

How to generate a random double number in the inclusive [0,1] range?

After taking a shower, I have conceived of a potential solution based on my understanding of how a random floating point generator works. My solution makes three assumptions, which I believe to be reasonable, however I can not verify if these assumptions are correct or not. Because of this, the following code is purely academic in nature, and I would not recommend its use in practice. The assumptions are as follows:

The distribution of random.NextDouble() is uniform
The difference between any two adjacent numbers in the range produced by random.NextDouble() is a constant epsilon e
The maximum value generated by random.NextDouble() is equal to 1 - e

Provided that those three assumptions are correct, the following code generates random doubles in the range [0, 1].

// For the sake of brevity, we'll omit the finer details of reusing a single instance of Random
var random = new Random();

double RandomDoubleInclusive() {
    double d = 0.0;
    int i = 0;

    do {
        d = random.NextDouble();
        i = random.Next(2);
    } while (i == 1 && d > 0)
    
    return d + i;
}

This is somewhat difficult to conceptualize, but the essence is somewhat like the below coin-flipping explanation, except instead of a starting value of 0.5, you start at 1, and if at any point the sum exceeds 1, you restart the entire process.

From an engineering standpoint, this code is a blatant pessimization with little practical advantage. However, mathematically, provided that the original assumptions are correct, the result will be as mathematically sound as the original implementation.

Below is the original commentary on the nature of random floating point values and how they're generated.

Original Reply:

Your question carries with it a single critical erroneous assumption: Your use of the word "Correct". We are working with floating point numbers. We abandoned correctness long ago.

What follows is my crude understanding of how a random number generator produces a random floating point value.

You have a coin, a sum starting at zero, and a value starting at one half (0.5).

Flip the coin.
If heads, add the value to the sum.
Half the value.
Repeat 23 times.

You have just generated a random number. Here are some properties of the number (for reference, 2^23 is 8,388,608, and 2^(-23) is the inverse of that, or approximately 0.0000001192):

The number is one of 2^23 possible values
The lowest value is 0
The highest value is 1 - 2^(-23);
The smallest difference between any two potential values is 2^(-23)
The values are evenly distributed across the range of potential values
The odds of getting any one value are completely uniform across the range
Those last two points are true regardless of how many times you flip the coin
The process for generating the number was really really easy

That last point is the kicker. It means if you can generate raw entropy (i.e. perfectly uniform random bits), you can generate an arbitrarily precise number in a very useful range with complete uniformity. Those are fantastic properties to have. The only caveat is that it doesn't generate the number 1.

The reason that caveat is seen as acceptable is because every other aspect of the generation is so damned good. If you're trying to get a high precision random value between 0 and 1, chances are you don't actually care about landing on 1 any more than you care about landing on 0.38719, or any other random number in that range.

While there are methods for getting 1 included in your range (which others have stated already), they're all going to cost you in either speed or uniformity. I'm just here to tell you that it might not actually be worth the tradeoff.

Use rand() to generate uniformly distributed floating point numbers on (a,b), [a,b), (a,b], and [a,b]

If you want every double in the range to be possible, with probability proportional to the difference between it and its adjacent double values, then it's actually really hard.

Consider the range [0, 1000]. There are an absolute bucketload of values in the very tiny first part of the range: a million of them between 0 and 1000000*DBL_MIN, and DBL_MIN is about 2 * 10^-308. There are more than 2^32 values in the range altogether, so clearly one call to rand() isn't enough to generate them all. What you'd need to do is generate the mantissa of your double uniformly, and select an exponent with an exponential distribution, and then fudge things a bit to ensure the result is in range.

If you don't require every double in the range to be possible, then the difference between open and closed ranges is fairly irrelevant, because in a "true" continuous uniform random distribution, the probability of any exact value occurring is 0 anyway. So you might as well just generate a number in the open range.

All that said: yes, your proposed implementations generate values that are in the ranges you say, and for the closed and half-closed ranges they generate the end-points with probability 1/(RAND_MAX+1) or so. That's good enough for many or most practical purposes.

Your fiddling around with +1 and +2 works provided that RAND_MAX+2 is within the range that double can represent exactly. This is true for IEEE double precision and 32 bit int, but it's not actually guaranteed by the C standard.

(I'm ignoring your use of long double because it confuses things a bit. It's guaranteed to be at least as big as double, but there are common implementations in which it's exactly the same as double, so the long doesn't add anything except uncertainty).