On Sat, Dec 21, 2019 at 11:29:20PM +0000, Alexander Nasonov wrote:
> 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.
A very clear way to generate an equidistributed floating point number in
the range [0.0, 1.0) given a number n in the range [0, 2^48-1] is simply:
return n * 0x1.0p-48
Requiring only C99 and <stdint.h>.
>
> --
> Alex
--
/ Raimo Niskanen, Erlang/OTP, Ericsson AB
