Random doubles distributed between -DBL_MAX and DBL_MAX

Question

I am trying to generate random doubles with a roughly equal probability for all representable finite values. I believe this would be something like a randomly signed exponential_distribution, a uniform_distribution won't work because the representable doubles are not uniformly distributed. I have this hacky code that seems to do what I want:

#include <cmath>
#include <iostream>
#include <random>

template < class GEN >
double double_distribution (GEN& generator)
{
    using gen_type = typename GEN::result_type;
    static_assert (sizeof (gen_type) == sizeof (double), "");
    union {gen_type g; double f;} result;
    do {
        result.g = generator();
    } while (!std::isfinite (result.f));
    return result.f;
}

int main() {
    std::mt19937_64 mt ((std::random_device())());
    for (int n = 0; n < 10; ++n) {
        std::cout << double_distribution (mt) << ' ';
    }
    std::cout << '\n';
}

but I'm wondering if there is a simpler way to approximate this using one of the STL's existing RandomNumberDistribution classes?

Answer 1

If you assume a 64-bit IEEE double, here's a way. (For a different representation of double, you should be able to modify this basic approach.)

The 64 bits are divided into 1 bit for sign, 11 bits for exponent, and 52 bits for the fraction. Special values like NaN and inf are represented by special values of the exponent bits. You have two basic choices from here:

1. Generate a random number from the set of all possible 64 bit combinations (like a uniform on an appropriately sized integer type), and set the bits on your double from that. In this case, you'll sometimes draw a number that corresponds to, e.g., NaN. You'll either need to discard those and redraw or accept them, depending on whether you want NaN returned. (Sounds like you probably don't in your case.)

2. Generate the bits for the sign and fractional part in one draw and then separately draw the exponent. If you do this, you should be able to modify your draw such that you don't pull NaN values - The special value seem to correspond to extreme values of the exponent, so this would just amount to changing the range of your draw on the integer, which I think will be supported easily by many packages. In this case, you don't have to "re-draw" to avoid the special values, but you probably end up with some "extra" bits on each draw since you'll probably end up pulling from a distribution for an integer type (with a multiple-of-8 number of bits), but using 53 for the sign+fraction and only 11 for the exponent. You could probably do something smart to economize these bits so that "thrown out" bits from one draw are used in the next draw. You'd have to think through to make sure that this doesn't introduce a dependence between your draws, but I think it should be ok.

If efficiency isn't a big deal for you, then #1 looks simpler to implement. If you need efficiency, then maybe #2 could be better, although probably at the expense of more complicated code.

Answer 2

What you're looking for is std::uniform_real_distribution.

Note that this will - at least on Visual Studio - not work for the entire range from numeric_limits<double>::lowest() to numeric_limits<double>::max() because the calculation involves the distance between the lower and upper bounds of the range, which is obviously larger than double allows for.

I'd suggest you to rethink your approach anyway, because it's probably not actually what you want - a uniform distribution between all double values is almost exclusively going to generate values with an extremely large magnitude (almost all numbers you'll get are going to be something-something*10^300).

In any case, if you actually want to use a uniform distribution here's an example:

#include <random>
#include <limits>
#include <iostream>
#include <cmath>

int main()
{
    std::random_device rd;
    std::mt19937 gen(rd());
    // sqrt so distance between min and max fits into a double
    auto max = std::sqrt(std::numeric_limits<double>::max());
    auto min = -max;
    std::cout << "min: " << min << " max: " << max << std::endl;
    std::uniform_real_distribution<double> dis(min, max);
    for (int n = 0; n < 100; ++n) {
        std::cout << dis(gen) << '\n';
    }
    std::cout << std::endl;
}