5 ms·

Fixing Random

44 points by azhenley 8y ago

azhenley 8y ago
Part 2: https://ericlippert.com/2019/02/04/fixing-random-part-2/ https://ericlippert.com/2019/02/04/fixing-random-part-2/ Part 3: https://ericlippert.com/2019/02/07/fixing-random-part-3/ https://ericlippert.com/2019/02/07/fixing-random-part-3/
deleted 8y ago
[deleted]
a1369209993 8y ago
> An alternative, and possibly better solution, would be to build the double from bits directly, rather than doing an expensive division. Clever trick: uint64 x = 0x3FF0'0000'0000'0000 | (rand64()&0xF'FFFF'FFFF'FFFF) return *(double*)&x - 1.0 This gives a perfect uniform distribution on [0,1), with the caveat that everything is rounded to a multiple of 2^-52 (≈ 2.22e-16).
- pjscott 8y ago
  I love this trick. It generates a float of the form 1.xxxxx..., filling in the x part with random bits, generating a random number between 1 and 2. Then it subtracts 1.
  - FabHK 8y ago
    Note though that the least significant bit of your float will always be zero, if I understand correctly (and the subsequent ones will be zero more often than not). PRNG are tricky, and floats are tricky, and the combination is tricky, unsurprisingly. See this discussion here: https://github.com/JuliaLang/julia/issues/16344 https://github.com/JuliaLang/julia/issues/16344
    - pjscott 8y ago
      Aw, man. Some day we'll have this all figured out, but that day won't be soon.
- kbsletten 8y ago
  So, am I correct that this would give you effectively 52 random bits instead of the maximum 53? Is the problem that the 53rd bit would mean there are possibly invalid values you'd have to throw out, or that the distribution wouldn't be uniform?
  - a1369209993 8y ago
    The problem is that: a, there are over a thousand random bits to fill, and you'll need a significantly more sophisticated mechanism to fill them; b, if you need more than a dozen significant digits, 64-bits floats probably aren't what you want anyway; and c, clever tricks aren't the place for rigorous correctness when "correct to within rounding errors" will do just as well. If you really want you could tack on: ret += (entropy >> 63) * 0x1.0p-53; // 1.110223024625156540e-16; OR ret += (entropy >> 52) * 0x1.0p-64; // 5.421010862427522170e-20; right before the return to pack in some extra bits.
    - kbsletten 8y ago
      This reply made me even more confused. According to my naive understanding of floating point there are 53 bits for the mantessa (sp?) so in my head, you would specify a exponent of 2^(-53) and fill the mantessa with random bits to get a "pretty good" random number. If you wanted to represent every possible value of a double, then you have to vary the exponent and your distribution becomes non-uniform. I assume I'm missing something fairly obvious to others, but your value seems to be exactly one bit short of the "ideal" (if you call it that, but really it's just my imagination of a "perfect" random double). I see your reply mentions you can get back that bit, but I don't understand why it has to be it's own step. What limitation of doubles are you overcoming when you use the range 1-2 and then subtract one instead of just generating in the range 0-1?
      - jacobolus 8y ago
        For IEEE floats, the range [1, 2) has uniform absolute precision. Not so for [0, 1), which spans many values of the exponent. There is uniform absolute precision on each interval among [0.5, 1), [0.25 0.5), [0.125, 0.25), &c., but as you get closer and closer to 0, there are still just as many fraction bits available (same relative precision), so these numbers are much more precise in absolute terms (i.e. relative to 1.0). You can generate a random number in [1, 2) by just filling in the fraction part with random bits, because the exponent is always the same. Then you can subtract 1 to get a number in the range [0, 1). But for a number in e.g. the range [0.125, 0.25), the last 3 bits of the fraction will always be 0, meaning you are skipping 7/8 of the possible floating point numbers in that range. The smallest possible denormal number among IEEE doubles is 2^(–1074), but the smallest non-zero value you will get with this method is 2^(–52). If you want to generate any possible floating point number in the range [0,1) with probability proportional to the difference between it and the next representable number, it gets pretty tricky / computationally expensive.
- deleted 8y ago
  [deleted]
- redcalx 8y ago
  See NextDoubleHighRes() here: https://github.com/colgreen/Redzen/blob/master/Redzen/Random/RandomSourceBase.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... There's an explanation in the comments, and a link to the original fuller explanation and source, here: https://mumble.net/~campbell/tmp/random_real.c https://mumble.net/~campbell/tmp/random_real.c I.e. I'm pretty sure your approach gives a non uniform distribution, but happy to discuss in more detail if you want.
  - FabHK 8y ago
    Among the number it emits, the distribution is "uniform", but it leaves out many possible doubles (half between 0.5 and 1, a quarter between 0.25 and 0.5, etc.), as far as I understand. Great sources! (EDIT: replace "omit" by "leave out", because it's too close to "emit")
    - redcalx 8y ago
      Yeh that sounds right. Had forgotten the details as it's been a while since I went through all the issues and wrote and tested the redzen code.
  - a1369209993 8y ago
    > your approach gives a non uniform distribution Unless I made some significant error, it is a uniform distribution over multiples of 2^-52. It's also a discrete approximation of the real uniform distribution in that for any subdivision of [0,1) into N equal (half-open) intervals, the most-populated interval has at most one more inhabitant the least populated one. It is however strictly inferior to rand53()*2^-53 (NextDoubleInner at your github link) unless your int64->double or float-multiply operations are quite bad.
iheartpotatoes 8y ago
I'm not following his reason for writing this. In part 2 he says he doesn't care about the performance hit from a crypto library that implements a system TRNG, sooooo, use that? https://developer.mozilla.org/en-US/docs/Web/API/Crypto/getRandomValues https://developer.mozilla.org/en-US/docs/Web/API/Crypto/getR... Or even study the thread on entropy pool paranoia from the Stanford Javascript Crypto Library: https://github.com/bitwiseshiftleft/sjcl/issues/77 https://github.com/bitwiseshiftleft/sjcl/issues/77 /scratches head at FAIC/
- detaro 8y ago
  Why are you linking to Javascript docs under an article about why the default C# APIs for random are bad? And even if there's already replacement libraries, a discussion of how they can be better than the default is useful?
- ericlippert 8y ago
  The reason for writing this is to gently introduce the general topic of improving how we represent stochastic workflows in mainstream programming languages. Entering that complex topic by exploring just how horrid even the most basic aspects of randomness are represented is an accessible and compelling entry point. This is part three of a very long series that will explore the connections between LINQ, probability distributions, and the difficulties of sampling from arbitrary distributions.
redcalx 8y ago
There are a number of issues with System.Random, all of the known issues to date are described in this open github issue: https://github.com/dotnet/corefx/issues/23298 https://github.com/dotnet/corefx/issues/23298 The redzen library contains RNGs that address all of the known issues, e.g. see: https://github.com/colgreen/Redzen/blob/master/Redzen/Random/RandomSourceBase.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... https://github.com/colgreen/Redzen/blob/master/Redzen/Random/Xoshiro256StarStarRandom.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... That Xoshiro256StarStarRandom source has also been merged into Math.Net.
SOLAR_FIELDS 8y ago
I personally fell victim to this very mistake he mentions when I first started with C#. It really doesn’t follow the principle of least astonishment and it makes me feel a bit better knowing that it’s apparently a pretty common mistake.

44 points by azhenley 8y ago

Part 2: https://ericlippert.com/2019/02/04/fixing-random-part-2/ https://ericlippert.com/2019/02/04/fixing-random-part-2/ Part 3: https://ericlippert.com/2019/02/07/fixing-random-part-3/ https://ericlippert.com/2019/02/07/fixing-random-part-3/

deleted 8y ago

[deleted]

a1369209993 8y ago

> An alternative, and possibly better solution, would be to build the double from bits directly, rather than doing an expensive division. Clever trick: uint64 x = 0x3FF0'0000'0000'0000 | (rand64()&0xF'FFFF'FFFF'FFFF) return *(double*)&x - 1.0 This gives a perfect uniform distribution on [0,1), with the caveat that everything is rounded to a multiple of 2^-52 (≈ 2.22e-16).

pjscott 8y ago
I love this trick. It generates a float of the form 1.xxxxx..., filling in the x part with random bits, generating a random number between 1 and 2. Then it subtracts 1.
- FabHK 8y ago
  Note though that the least significant bit of your float will always be zero, if I understand correctly (and the subsequent ones will be zero more often than not). PRNG are tricky, and floats are tricky, and the combination is tricky, unsurprisingly. See this discussion here: https://github.com/JuliaLang/julia/issues/16344 https://github.com/JuliaLang/julia/issues/16344
  - pjscott 8y ago
    Aw, man. Some day we'll have this all figured out, but that day won't be soon.
kbsletten 8y ago
So, am I correct that this would give you effectively 52 random bits instead of the maximum 53? Is the problem that the 53rd bit would mean there are possibly invalid values you'd have to throw out, or that the distribution wouldn't be uniform?
- a1369209993 8y ago
  The problem is that: a, there are over a thousand random bits to fill, and you'll need a significantly more sophisticated mechanism to fill them; b, if you need more than a dozen significant digits, 64-bits floats probably aren't what you want anyway; and c, clever tricks aren't the place for rigorous correctness when "correct to within rounding errors" will do just as well. If you really want you could tack on: ret += (entropy >> 63) * 0x1.0p-53; // 1.110223024625156540e-16; OR ret += (entropy >> 52) * 0x1.0p-64; // 5.421010862427522170e-20; right before the return to pack in some extra bits.
  - kbsletten 8y ago
    This reply made me even more confused. According to my naive understanding of floating point there are 53 bits for the mantessa (sp?) so in my head, you would specify a exponent of 2^(-53) and fill the mantessa with random bits to get a "pretty good" random number. If you wanted to represent every possible value of a double, then you have to vary the exponent and your distribution becomes non-uniform. I assume I'm missing something fairly obvious to others, but your value seems to be exactly one bit short of the "ideal" (if you call it that, but really it's just my imagination of a "perfect" random double). I see your reply mentions you can get back that bit, but I don't understand why it has to be it's own step. What limitation of doubles are you overcoming when you use the range 1-2 and then subtract one instead of just generating in the range 0-1?
    - jacobolus 8y ago
      For IEEE floats, the range [1, 2) has uniform absolute precision. Not so for [0, 1), which spans many values of the exponent. There is uniform absolute precision on each interval among [0.5, 1), [0.25 0.5), [0.125, 0.25), &c., but as you get closer and closer to 0, there are still just as many fraction bits available (same relative precision), so these numbers are much more precise in absolute terms (i.e. relative to 1.0). You can generate a random number in [1, 2) by just filling in the fraction part with random bits, because the exponent is always the same. Then you can subtract 1 to get a number in the range [0, 1). But for a number in e.g. the range [0.125, 0.25), the last 3 bits of the fraction will always be 0, meaning you are skipping 7/8 of the possible floating point numbers in that range. The smallest possible denormal number among IEEE doubles is 2^(–1074), but the smallest non-zero value you will get with this method is 2^(–52). If you want to generate any possible floating point number in the range [0,1) with probability proportional to the difference between it and the next representable number, it gets pretty tricky / computationally expensive.
deleted 8y ago
[deleted]
redcalx 8y ago
See NextDoubleHighRes() here: https://github.com/colgreen/Redzen/blob/master/Redzen/Random/RandomSourceBase.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... There's an explanation in the comments, and a link to the original fuller explanation and source, here: https://mumble.net/~campbell/tmp/random_real.c https://mumble.net/~campbell/tmp/random_real.c I.e. I'm pretty sure your approach gives a non uniform distribution, but happy to discuss in more detail if you want.
- FabHK 8y ago
  Among the number it emits, the distribution is "uniform", but it leaves out many possible doubles (half between 0.5 and 1, a quarter between 0.25 and 0.5, etc.), as far as I understand. Great sources! (EDIT: replace "omit" by "leave out", because it's too close to "emit")
  - redcalx 8y ago
    Yeh that sounds right. Had forgotten the details as it's been a while since I went through all the issues and wrote and tested the redzen code.
- a1369209993 8y ago
  > your approach gives a non uniform distribution Unless I made some significant error, it is a uniform distribution over multiples of 2^-52. It's also a discrete approximation of the real uniform distribution in that for any subdivision of [0,1) into N equal (half-open) intervals, the most-populated interval has at most one more inhabitant the least populated one. It is however strictly inferior to rand53()*2^-53 (NextDoubleInner at your github link) unless your int64->double or float-multiply operations are quite bad.

iheartpotatoes 8y ago

I'm not following his reason for writing this. In part 2 he says he doesn't care about the performance hit from a crypto library that implements a system TRNG, sooooo, use that? https://developer.mozilla.org/en-US/docs/Web/API/Crypto/getRandomValues https://developer.mozilla.org/en-US/docs/Web/API/Crypto/getR... Or even study the thread on entropy pool paranoia from the Stanford Javascript Crypto Library: https://github.com/bitwiseshiftleft/sjcl/issues/77 https://github.com/bitwiseshiftleft/sjcl/issues/77 /scratches head at FAIC/

detaro 8y ago
Why are you linking to Javascript docs under an article about why the default C# APIs for random are bad? And even if there's already replacement libraries, a discussion of how they can be better than the default is useful?
ericlippert 8y ago
The reason for writing this is to gently introduce the general topic of improving how we represent stochastic workflows in mainstream programming languages. Entering that complex topic by exploring just how horrid even the most basic aspects of randomness are represented is an accessible and compelling entry point. This is part three of a very long series that will explore the connections between LINQ, probability distributions, and the difficulties of sampling from arbitrary distributions.

redcalx 8y ago

There are a number of issues with System.Random, all of the known issues to date are described in this open github issue: https://github.com/dotnet/corefx/issues/23298 https://github.com/dotnet/corefx/issues/23298 The redzen library contains RNGs that address all of the known issues, e.g. see: https://github.com/colgreen/Redzen/blob/master/Redzen/Random/RandomSourceBase.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... https://github.com/colgreen/Redzen/blob/master/Redzen/Random/Xoshiro256StarStarRandom.cs https://github.com/colgreen/Redzen/blob/master/Redzen/Random... That Xoshiro256StarStarRandom source has also been merged into Math.Net.

SOLAR_FIELDS 8y ago

I personally fell victim to this very mistake he mentions when I first started with C#. It really doesn’t follow the principle of least astonishment and it makes me feel a bit better knowing that it’s apparently a pretty common mistake.