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

[sunoru]

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