| Subject: | rand is only filling 32 bits |
The rand() function only returns 32 bits of random data. This is done in _xorshift.c:
return upper_limit * ((double)xorshift_irand(prng) / UINT32_MAX);
While this is fast, it gives lousy results. We should be filling 53 bits. With 64-bit types we could do:
return (irand64() >> 11) * (1.0 / (1UL << 53));
or stick with 32-bit and do something like:
uint32_t a = xorshift_irand(pong) >> 6;
uint32_t b = xorshift_irand(pong) >> 5;
return upper_limit * (a * 134217728.0 + b) / 9007199254740992.0;
or the method used by the MTwist module which creates and caches the constants to use.
Note the interval chosen by the API needs to be considered as well. Typically rand() uses the half-open interval [0,1), while the function currently implemented uses [0,1] as far as I can tell.