This was pretty much my reasoning as well, but I missed that you can generate a floating-point integer in [0,2^53) and then scale down to [0,1) with a multiply,

This was pretty much my reasoning as well, but I missed that you can generate a floating-point integer in [0,2^53) and then scale down to [0,1) with a multiply, which gets you one additional significant bit for free or cheaper, and simplifys the code. Edit: I think for the general case you want something like: uint64 x = rand64() & 0x001F'FFFF'FFFF'FFFF double xf = .5 * 1p-52 * x // return xf; /* 53 bits of entropy */ if(xf < .5) { y = rand64() & 0x000F'FFFF'FFFF'FFFF double yf = .5 * 1p-52 * 1p-52 * y xf += yf if(xf < .5 * 1p-52) ... /* loop about 20 times */ }