art4711 / random-double Goto Github PK

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

(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.