art4711 / random-double

An algorithm for generating random doubles.

Generating random doubles

The first step was to figure out how to generate numbers in the range [0.0,1.0). Those experiments and the reasoning behind them is in comments and code in rd.c.

A few months later I looked at gcc and llvm standard c++ libraries to see if their std::uniform_real_distribution actually solved this problem correctly. They haven't. The result of that is in urd.cxx. The bug report for llvm is here, I haven't made a bug report for gcc.

This triggered me to actually figure out how to extend this to arbitrary ranges. The first attempt is documented in comments and code in arbitrary_range.c, but it only deals with positive numbers for now. The nice thing about it is that despite a completely different algorithm the generated numbers from arbitrary_range.c set the exact same bits as the numbers in rd.c which means that either my early assumptions were correct or that my brain farts are at least consistent.

TODO

• Tackle negative numbers. Naively it should just be like arbitrary_range, except that we find the biggest absolute value of (from,to), use that to find the step and count and do some sign flipping in the right place. But that's for another day.

DISCLAIMER

Please notice that I'm not claiming that anything above makes sense or that the tests are even close to being comprehensive. A programmer who tests his own edge cases is like an author writing his own book reviews.

To be useful this would require vastly more testing and independent verificaiton from someone who actually understands set theory and ieee 754 floating point numbers. I can confidently say that I know neither, but as the industry standard on stackoverflow and the c++ libraries seems to indicate, neither do most other people.

random-double's Issues

The `r0to1` functions actually don't return `0.0` with a 2^-53 probability.

double
r0to1(void)
{
    int e;
    uint64_t m;

    e = ffsll(rX(52));
    if (e == 0)
        return 0.0;
    m = rX(52 - e + 1) << (e - 1);
    return ldexp(0x1p52 + m, -52 - e);
}

double
r0to1b(void)
{
    uint64_t r = rX(53);
    int e = ffsll(r);
    uint64_t m;
    if (e > 52 || e == 0)
        return 0.0;
    /* Shift out the bit we don't want set. */
    m = (r >> e) << (e - 1);
    return ldexp(0x1p52 + m, -52 - e);
}

The probability should be 2^-52 for a 52-bit uint x to have ffsll(x) return 0, and 2^-53 for a 53-bit uint x to have ffsll(x) return 53 or 0 (which add up to 2^-52).

The functions can equiprobably (2^-53) generate all the possible return values, except for 0 and 2^-53. The latter one isn't returned, making it more possible for 0. In fact just removing the e > 52 condition in r0to1b will make it:

double
r0to1c(void)
{
    uint64_t r = rX(53);
    int e = ffsll(r);
    uint64_t m;
    if (e == 0)
        return 0.0;
    m = (r >> e) << (e - 1);
    return ldexp(0x1p52 + m, -52 - e);
}

In this case, when e == 53, m will be 0, and then the return value will be ldexp(0x1p52, -105), i.e. 2^-53.

btw I have another question: have you ever done some benchmark for the different approaches to convert uint64_t to double? TIA

Relevant prior art: dSFMT

(I tried to post this as a comment to your blog post, but it didn't post, and there was no message telling me it was going into a moderation queue. So I apologize if this is redundant.)

You might want to take a look at dSFMT, which directly generates double-precision random numbers, including uniform on [0, 1): http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/

This algorithm still has undesirable statistical properties, however, as discussed at JuliaLang/julia#6464.