Ingo Schwarze wrote:
> Looking at our code in lib/libc/stdlib/drand48.c, i conclude that
> drand48(3) does return 0.0 with a probability of 2^-48.
I looked at the code too and I have some comments.
> More generally, the function returns a uniform distribution of
> numbers from the set {2^-48 * n | n integer and 0 <= n < 2^48}.
You don't need three ldexp calls to compose 2^-48 * n:
uint64_t n = (uint64_t)xseed[2] << 32 | xseed[1] << 16 | xseed[0];
return ldexp((double)n, -48);
> Talking about loading bits into the mantissa and adjusting the
> exponent feels mildly confusing, given that the distribution is
> simply uniform over a fixed finite set. I'm not sending a patch
> because i never looked at how floating point representation works
> internally, so i would likely only make matters worse.
Not sure if you wanted to patch the man page or the code.
If the latter, take a binary representation of 1.0 and replace 52
least significant bits with (pseudo)random bits then subtract 1.0.
double d = 1.0;
uint64_t u = 0;
memcpy(&u, &d, sizeof(d));
u |= n << 4; /* scale 48 bits to 52 bits, n is from the previous
snippet */
memcpy(&d, &u, sizeof(u));
return d - 1.0;
It will only work with IEEE 754 compliant doubles, though.
--
Alex
