Replace the rand() with a portable rng

[aegkmq]

close #123 , close #21

[xefoci7612]

Sorry, I don't understand this line:

return (random_u32() >> 9) * (1.0f / 8388608.0f);

The conventional/idiomatic way to get a random float from a random integer seems to be:

float x = (float)rand()/(float)RAND_MAX_64

where in our case (64 bit): RAND_MAX_64 = 0x3FFFFFFF

Maybe we can just use the original algorithm (without the last artificial right shift), and keep the conversion to float explicit as is now. Something like:

#define RAND_MAX_64  0x3FFFFFFFULL

unsigned long long rng_seed;
unsigned int random_u64() {
    rng_seed ^= rng_seed >> 12;
    rng_seed ^= rng_seed << 25;
    rng_seed ^= rng_seed >> 27;
    return rng_seed * 0x2545F4914F6CDD1Dull;
}
.....

    // sample index from probabilities, they must sum to 1
    float r = (float)random_u64() / (float)RAND_MAX_64;

[aegkmq]

This paper says that to obtain a float in [0, 1), we should generate a random integer in the interval [0, 2^p - 1] and divide it by 2^p. (p=53 for double, p=24 for float)
So I'll modify p and keep the current method.

[karpathy]

tbh it's a little bit scary to me that we're just winging it in a random PR. Makes me uncomfortable about correctness. Are there any verifiable references to some reputable repos that use this code nearly verbatim?

[aegkmq]

I copied that conversion code from a previous version of numpy. (They corrected the divisor for the same reason)

[karpathy]

Got it, thank you @aegkmq . I like this, I think I'll merge it later this evening.

[m1lkweed]

The conventional/idiomatic way to get a random float from a random integer seems to be:

float x = (float)rand()/(float)RAND_MAX_64

As #123 mentions, the idiomatic way isn't really good for generating random floats. I wrote this up in godbolt as an example. Numbers close to 1 are rounded to 1 and numbers near 0 drop down near 2e-10 with the idiomatic code. rand_float() produces more sensible values.

[m1lkweed]

And here is an example with the code from this issue (return (a >> 8) / 16777216.0f;)

[aegkmq]

It's just the difference between truncation methods. Your example increases i by one, so & 0x7fffff has an unfair advantage on small i versus >> 8. rand_float returns the duplicate and non-decreasing values on larger i. It only indicates a difference between truncation methods and does not mean one method is wrong. Also, (a & 0xffffff) / 16777216.0f produces more sensible results than rand_float in your test.

[karpathy]

Thank you all, I'll take it.

[cshenton]

Just FYI, xoshiro/xoroshiro PRNGs (such as the one I PR'd earlier in the week) are both faster and statistically more sound than xorshift PRNGs. Details: https://prng.di.unimi.it/ and https://xoshiro.di.unimi.it/hwd.php

If for whatever reason you happen to prefer a full 64 bit PRNG, even if 40 of those bits will be thrown away such as they are here, I'd at least suggest using xoshiro**/xoshiro++ PRNGs, which are now the default PRNG in .NET 6, Lua, Zig, and GNU Fortran.

Though, in the interest of balance, Javascript does use a Xorshift PRNG.