That Wiki page makes it sound very complicated but for this purpose our implementation can be laughably simple which has the advantage that you know why it works and can maintain it properly with confidence.
Get suitably large inputs, for example if you're trying to pick integers between 2 and 11 inclusive, a nibble (half a byte) would be fine. Now, is the random input in the range you wanted? If so, you've got your answer. If not, throw this random input away and get more.
Too many programmers act as though random numbers were a precious resource.
Random numbers aren't precious, but they do take time to generate. Rejection sampling adds a branch and it can matter how often it's taken. In my property-testing framework, generating a large, random array of numbers more efficiently improved performance.
If your RNG outputs a u64 (a getrandom() call can do this, or even a simple RNG like xoshiro), a random int in TFA's range of [0, 10) has a very small (6 in 2⁶⁴ chance, or 3.2e-17%, might as well be 0) chance of needing to redraw.
Obviously, wider ranges will (possibly) reject more often, but unless the range is huge, rejection should be a very cold branch.
You’re throwing a lot of bits away, though. If you used four bits per number then you could pack 16 random digits into one 64-bit draw. But a rejection-sampling branch will be taken much more often.
If you want to split a 64-bit random number into 16 4-bit random numbers, the quality of the RNG that generates the 64-bit random numbers must be extremely high. Only cryptographic RNGs may have such a high quality.
Any non-cryptographic PRNG is useful only if it is much faster than a cryptographic RNG, e.g. one using AES, which in many modern CPUs needs around 10 clock cycles to generate a 128-bit random number (less than that, even less than half of that, in some more recent CPUs).
So non-cryptographic PRNGs must generate a 64-bit number in less than 1 nanosecond to be competitive. Many older PRNGs are not this fast, so they are completely obsolete.
Such PRNGs do not offer any guarantee that if you take more than one piece of a generated number they will not be correlated. So the rule is that you can not make 2 or more random numbers from 1 random number provided by a PRNG.
Rejection is unavoidable if you want a uniform distribution for a number of choices that is not a power of two.
Nonetheless, the rejection can be done either before or after the multiplication or division that does most of the job for passing from the input range to the output range.
The place where rejection is done can be chosen to minimize the amount of values that are rejected, making thus unlikely that the branch is taken, so it will be correctly predicted most of the time.
For maximum speed, it is preferable to use multiplication instead of division, i.e. the input is seen as a fraction less than 1 and after multiplication only the integer part is retained. Taking care to use multiplication is normally more important than worrying that rejection may decrease the performance.
Wow that page really gets lost in the weeds of multiple dimensions.
And yeah it's just rerolling when your random number is out of range. If you want it as simple as possible, always generate from 0-n, and grab barely enough random bits for n to fit.
That's what for example the Doom "fizzle" effect does. You can pick taps for an LFSR that generates a sequence 2^n -1 steps long - it can never contain zero and all numbers from 1 to 2^n -1 will appear exactly once.
So to "fizzle" a 320x200 resolution you generate a 17-bit LFSR, the lower 8 bits being the Y co-ordinate and the upper 9 being the X co-ordinate (or the other way round, don't you let me tell you how to live). Some of those are off the screen! Ignore them. The "off the edges" fizzles that you don't draw are so randomly dispersed that you never actually notice them.
Then eventually the LFSR loops back to the seed value, and you know you've done them all.
You can always choose not to reject again if some number of rejections fail. With a sufficiently large random number range, the chances of rejection at any step can be made tiny, and the possibility of successive rejections drops exponentially with the number of rejections allowed until you have a theoretical probability smaller than the probability of hardware bit errors on the platform itself. If this happens, you can either throw once more and use that result, with a tiny amount of modulo bias, or you can just throw an exception.
That chance is strictly zero though. While it's theoretically possible for the loop to never terminate, such an event is basically never going to happen in any real world scenario.
If you reduce the range to the nearest power of two then it's expected to only need 2 iterations on average.
And if the number generator is guaranteed to generate every bit pattern eventually, then it's guaranteed to terminate.
Using 10 integers was actually a simplification to explain what you're talking about here.
If the basic available function was rand_1_to_10(), modulo-ing its result by 3 would clearly make one result more probable than the others. The same happens with random_u64(); even though the difference is smaller, it can still be significant.
Assuming the code is black box, an outside observer would need no more than a thousand observations before they could definitely say something was off about your rand_1_to_10() mod 3 code (with 95% confidence).
But in the case of the random_u64() mod 10, you would need the heat death of the universe number of samples before that bias would be noticeably distinguishable from random.
Is there an actual example of where that specific case could be significant for anything other than maybe nation-state cryptography?
this is what has always baffled me about statistics...having the opportunity to make things exact, with little effort, it is dismissed just because the small error
That sounds related to the pigeon hole principle: if you have n objects and m buckets, and m doesn’t divide n, then some buckets have more objects than others.
A Uniform(0, 1) PRNG is the right primitive on which to base a PRNG library for arbitrary real-valued random variables because any real-valued random variable can be represented as an inverse quantile transform of a Uniform(0, 1) random variable. While I agree that only providing this primitive risks footguns as described in the article, I'm skeptical that providing a better UI alone would meaningfully reduce the risk of such footguns because availability and convenience is no guarantee of use if a prevailing attitude of users is that they know better and don't need it. Bisection search is simpler than rejection sampling or inverse transform sampling, yet it's common to see buggy, hand-rolled implementations of bisection search despite wide availability of library implementations with better UI ergonomics than PRNGs. I think wider use of fuzz testing or deterministic simulation testing really is necessary to disabuse people of that notion, along with more articles like the above explaining why hand-rolling an adapter to a uniform PRNG is a false economy compared to proven implementations with vetted statistical properties.
PRNG is the answer, if you want over a set of small numbers (say a 100) to get a nice distribution of your choice , then you shouldn’t use a RNG, you should generate the desired points on the curve and then shuffle them (or use deterministic code to that effect). A true RNG can output any numbers and only over a large set of numbers the distribution would be uniform.
> I'm skeptical that providing a better UI alone would meaningfully reduce the risk of such footguns because availability and convenience is no guarantee of use if a prevailing attitude of users is that they know better and don't need it.
If we transfer that argument to other domains it quickly falls apart: We don't need a safety on handguns since some people will opt to not use it properly. We don't need safety belts in cars since some people won't use it. We don't need handrails..
You get the point. This is basically the Nirvana fallacy (also called the perfect solution fallacy): a safety measure not being able to catch 100% of the cases it was meant to prevent is not an argument against it. We could have an argument if that measure would significantly impact usual use cases. Any safety measure needs to be judged both by the benefits and by how practical it is to deploy it (aside from other considerations like maintenance, etc)
In this case having an easier to use API is the opposite of impractical. It both safes developers time and reduces the number of errors. And if you still need to write your own function, you totally can. To me that sounds as close as you can get to the definition of a "no-brainer".
Yes, I should probably clarify and amend my argument.
My impression of the article is that it suggests if ergonomic APIs were more available than they already are, then the frequency with which we see bugs in implementing random variable sampling would decrease. I question that premise because the ergonomic APIs already widely exist -- as the article itself points out, in many language standard libraries, and also in popular third-party libraries. Such a suggestion seems like it relies too much on a single mitigation.
To borrow your examples:
- we don't rely on safety belts alone to reduce injuries: we augment that control with more engineering controls, regulations, and education.
- we don't rely on gun safeties alone to reduce injuries: we augment that control with more engineering controls, regulations, and education.
For obvious reasons, I think regulation or mandatory education would be impractical for this situation.
Because other engineering domains, including the ones you mentioned, rely on complementary controls to achieve further safety, I speculated that fuzz testing and deterministic simulation testing could be effective potential complementary controls. Property testing could be another.
I make this argument because I work as a computational statistician with software engineers in a large software company and I see variations on the theme of these bugs on a weekly basis. The usual justification I get is "I know that (insert library implementation) exists but I thought what I wrote was a simpler version of the same thing, and then I have fewer dependencies", as if equivalence of function were self-evident. It's not, and sometimes I can persuade the engineer to use the library from first principles, but generating witnesses violating the properties of the object they're computing tends to be more compelling, and generalizable strategies for testing tend to get more buy-in because the value-add is not limited to statistical applications.
There is a series of blog posts[0] by Eric Lippert called "Fixing Random"[1] that starts at random and covers distributions, weights, noise and so much more. It's a fantastic deep dive into the subject.
IIRC, Daniel Lemire has another better article with a more complete discussion of the methods for obtaining a uniform distribution when the number of choices is not a power of two.
The article linked by you is slightly misleading, because this paragraph is not really right:
> If you are willing to suffer a slight statistical bias, you can generate a floating-point values in [0,1) and multiply the result by the size of the interval. It avoids the division but introduces other overhead.
If you want to use multiplication instead of division there is no need whatsoever to generate random floating-point numbers.
You just interpret the random integer as a fixed-point fraction and you must use the unsigned integer multiply instruction that provides a double length product (or in some ISAs you must use the instruction that provides the high half of the product), from which you can extract the bits corresponding to the integer part of the result.
Thus there is no overhead in comparison with the division method. You should just use for input a PRNG for which all 2^N values are equiprobable, which is the case for most of them (there are some ancient PRNGs that provide only the non-negative signed integers; such PRNGs are wasting a bit and they must be shifted left by one before multiplication).
Moreover, there is absolutely no need to be "willing to suffer a slight statistical bias". If you apply rejection either before or after multiplication, you can remove any statistical bias, exactly like when you use division with rejection. You just need to be careful how you implement rejection, because the rejection criteria are different than for the division method. E.g., if you do the rejection after multiplication, you reject a result if its fractional part is greater than a constant that has been computed as the fractional part of the product between the multiplicand (i.e. the number of output values) and the greatest possible input value (2^N-1 interpreted as a fraction). This rejection ensures that the same number of input values fall in each output interval of unit length. When the number of output values is small in comparison with the number of input values, rejection happens very seldom.
While not saying this clearly, the "nearlydivisionless" function actually does exactly what I have said above, after reading just that introductory paragraph, which shows a straw man for comparison, but it is both more complex and slower than the multiplication method implemented correctly, which never requires any division. The "nearlydivisionless" function also removes the statistical bias by rejecting the result based on its fractional part, but instead of comparing with a precomputed constant it uses an inefficient method of computing the rejection threshold during function invocations, and using a division, even if the threshold can also be computed with a multiplication.
Moreover, I believe that the right implementation should use inline assembly for the multiply instruction, which provides guaranteed performance. The "nearlydivisionless" function avoids the inline assembly by writing some excessively complicated expression, with the hope that a clever compiler will optimize it by simplifying it to something equivalent with the correct assembly instruction. I do not believe that relying on very clever compilers is the right choice for achieving the portability of a program. On the contrary, even when recompiling the program for the same ISA, a new version of the same compiler may not do the desired optimization, causing unpredictable performance problems. If the compiler fails to optimize the 128-bit multiplication written in the "nearlydivisionless" function, by noticing that both 128-bit operands have been recently promoted from 64-bit values, it would do 4 multiplications instead of 1 multiplication.
Like I have said, I think that Daniel Lemire has another more recent article, including better multiplication methods, but I am too lazy to search for it now.
Wouldn't it be better multiplying the total of the choices with a random double. That way the random double is always strictly below total and the user isn't bitten by dumb float addition errors?
The only error then is when there are no choices or they set the weights to 0, negative or NaN. I am but a lowly musician though. I can't say I always think floats are simple...
I don’t think the “random uint” api is too low-level, or lacks a pit of success - I think you’re just reaching for the wrong api. The problem of “make n bits pseudo randomly set to either 1 or 0” is nearby to your problem of “choose an element according to a probability distribution”, but it’s a separate problem in its own right.
I think actually this is an argument for language designers to include a “std.choice” in their standard library that consumes random bytes and correctly performs common ergonomic operations like “get one element at random from this collection”.
(If your standard library tries to make a distinction between “regular random number generators” and “cryptographically secured random number generators”, I think this distinction between “generate random bits” and “make probabilistic choices” is about equally important.)
After finding out that trig/transcendental was basically not guaranteed to be equivalent across kernels (libc/musl), which caused the dreaded “only fails in CI” problem for me when I was trying to generate nonflat distributions of drng’s, I ended up creating https://github.com/pmarreck/random
to solve it, which it did
Ziggurat is badass. It’s also one of those algorithms (like “Metropolis-Hastings”) that has a fantastic name to go along with cool mathematical and practical properties.
I think this is interesting theory and fun to read about. If I was actually working with cryptography this would seem immensely important. I would also be arguing vehemently online about the std implementation of Mersenne twister and worrying about people analyzing bulk traffic with advanced scripts scraping a hundredth of a bit per sample.
My comment was mainly a joke. Since my work is primarily with performance and graphics, I rarely use any proper RNG and usually only need cheap things which look random enough at first glance.
Exactly. There's many a page on random which explains why taking the modulus of rand is the wrong choice, but it's also the easiest one, which makes it all to common.
A lot of crypto people bang on and on about choosing a proper random number generator, but if your inadequate random number generator is 1/10 less random than SOA, isn't RSA 2048 still unbreakable?
That Wiki page makes it sound very complicated but for this purpose our implementation can be laughably simple which has the advantage that you know why it works and can maintain it properly with confidence.
Get suitably large inputs, for example if you're trying to pick integers between 2 and 11 inclusive, a nibble (half a byte) would be fine. Now, is the random input in the range you wanted? If so, you've got your answer. If not, throw this random input away and get more.
Too many programmers act as though random numbers were a precious resource.
Obviously, wider ranges will (possibly) reject more often, but unless the range is huge, rejection should be a very cold branch.
Any non-cryptographic PRNG is useful only if it is much faster than a cryptographic RNG, e.g. one using AES, which in many modern CPUs needs around 10 clock cycles to generate a 128-bit random number (less than that, even less than half of that, in some more recent CPUs).
So non-cryptographic PRNGs must generate a 64-bit number in less than 1 nanosecond to be competitive. Many older PRNGs are not this fast, so they are completely obsolete.
Such PRNGs do not offer any guarantee that if you take more than one piece of a generated number they will not be correlated. So the rule is that you can not make 2 or more random numbers from 1 random number provided by a PRNG.
Nonetheless, the rejection can be done either before or after the multiplication or division that does most of the job for passing from the input range to the output range.
The place where rejection is done can be chosen to minimize the amount of values that are rejected, making thus unlikely that the branch is taken, so it will be correctly predicted most of the time.
For maximum speed, it is preferable to use multiplication instead of division, i.e. the input is seen as a fraction less than 1 and after multiplication only the integer part is retained. Taking care to use multiplication is normally more important than worrying that rejection may decrease the performance.
And yeah it's just rerolling when your random number is out of range. If you want it as simple as possible, always generate from 0-n, and grab barely enough random bits for n to fit.
So to "fizzle" a 320x200 resolution you generate a 17-bit LFSR, the lower 8 bits being the Y co-ordinate and the upper 9 being the X co-ordinate (or the other way round, don't you let me tell you how to live). Some of those are off the screen! Ignore them. The "off the edges" fizzles that you don't draw are so randomly dispersed that you never actually notice them.
Then eventually the LFSR loops back to the seed value, and you know you've done them all.
See also: https://romailler.ch/2020/07/28/crypto-modulo_bias_guide/
https://man.openbsd.org/arc4random.3#arc4random_uniform
random_u64() Mod 3 does indeed have a single bucket that is oversized. This overweights one option by about 5×10^-20.
rand() Itself has only 32767 possible values, so it's also common for a bucket to be overweighted depending on the number of buckets.
If the basic available function was rand_1_to_10(), modulo-ing its result by 3 would clearly make one result more probable than the others. The same happens with random_u64(); even though the difference is smaller, it can still be significant.
But in the case of the random_u64() mod 10, you would need the heat death of the universe number of samples before that bias would be noticeably distinguishable from random.
Is there an actual example of where that specific case could be significant for anything other than maybe nation-state cryptography?
For MingW maybe (due to MSVCRT). I think most libraries such as glibc have RAND_MAX at 2147483647.
https://en.wikipedia.org/wiki/Pigeonhole_principle
Most cases I have seen is a version of
Most people think "I understand this and I can use this" and just throw on top whatever they need beyond that. I was most people too.> I'm skeptical that providing a better UI alone would meaningfully reduce the risk of such footguns because availability and convenience is no guarantee of use if a prevailing attitude of users is that they know better and don't need it.
If we transfer that argument to other domains it quickly falls apart: We don't need a safety on handguns since some people will opt to not use it properly. We don't need safety belts in cars since some people won't use it. We don't need handrails..
You get the point. This is basically the Nirvana fallacy (also called the perfect solution fallacy): a safety measure not being able to catch 100% of the cases it was meant to prevent is not an argument against it. We could have an argument if that measure would significantly impact usual use cases. Any safety measure needs to be judged both by the benefits and by how practical it is to deploy it (aside from other considerations like maintenance, etc)
In this case having an easier to use API is the opposite of impractical. It both safes developers time and reduces the number of errors. And if you still need to write your own function, you totally can. To me that sounds as close as you can get to the definition of a "no-brainer".
My impression of the article is that it suggests if ergonomic APIs were more available than they already are, then the frequency with which we see bugs in implementing random variable sampling would decrease. I question that premise because the ergonomic APIs already widely exist -- as the article itself points out, in many language standard libraries, and also in popular third-party libraries. Such a suggestion seems like it relies too much on a single mitigation.
To borrow your examples:
- we don't rely on safety belts alone to reduce injuries: we augment that control with more engineering controls, regulations, and education. - we don't rely on gun safeties alone to reduce injuries: we augment that control with more engineering controls, regulations, and education.
For obvious reasons, I think regulation or mandatory education would be impractical for this situation.
Because other engineering domains, including the ones you mentioned, rely on complementary controls to achieve further safety, I speculated that fuzz testing and deterministic simulation testing could be effective potential complementary controls. Property testing could be another.
I make this argument because I work as a computational statistician with software engineers in a large software company and I see variations on the theme of these bugs on a weekly basis. The usual justification I get is "I know that (insert library implementation) exists but I thought what I wrote was a simpler version of the same thing, and then I have fewer dependencies", as if equivalence of function were self-evident. It's not, and sometimes I can persuade the engineer to use the library from first principles, but generating witnesses violating the properties of the object they're computing tends to be more compelling, and generalizable strategies for testing tend to get more buy-in because the value-add is not limited to statistical applications.
[0] https://ericlippert.com/2019/01/31/fixing-random-part-1/
[1] https://ericlippert.com/tag/fixing_random/
For sampling from a distribution, the best method is Vose's Alias. https://www.keithschwarz.com/darts-dice-coins/ is my favorite write up.
The article linked by you is slightly misleading, because this paragraph is not really right:
> If you are willing to suffer a slight statistical bias, you can generate a floating-point values in [0,1) and multiply the result by the size of the interval. It avoids the division but introduces other overhead.
If you want to use multiplication instead of division there is no need whatsoever to generate random floating-point numbers.
You just interpret the random integer as a fixed-point fraction and you must use the unsigned integer multiply instruction that provides a double length product (or in some ISAs you must use the instruction that provides the high half of the product), from which you can extract the bits corresponding to the integer part of the result.
Thus there is no overhead in comparison with the division method. You should just use for input a PRNG for which all 2^N values are equiprobable, which is the case for most of them (there are some ancient PRNGs that provide only the non-negative signed integers; such PRNGs are wasting a bit and they must be shifted left by one before multiplication).
Moreover, there is absolutely no need to be "willing to suffer a slight statistical bias". If you apply rejection either before or after multiplication, you can remove any statistical bias, exactly like when you use division with rejection. You just need to be careful how you implement rejection, because the rejection criteria are different than for the division method. E.g., if you do the rejection after multiplication, you reject a result if its fractional part is greater than a constant that has been computed as the fractional part of the product between the multiplicand (i.e. the number of output values) and the greatest possible input value (2^N-1 interpreted as a fraction). This rejection ensures that the same number of input values fall in each output interval of unit length. When the number of output values is small in comparison with the number of input values, rejection happens very seldom.
While not saying this clearly, the "nearlydivisionless" function actually does exactly what I have said above, after reading just that introductory paragraph, which shows a straw man for comparison, but it is both more complex and slower than the multiplication method implemented correctly, which never requires any division. The "nearlydivisionless" function also removes the statistical bias by rejecting the result based on its fractional part, but instead of comparing with a precomputed constant it uses an inefficient method of computing the rejection threshold during function invocations, and using a division, even if the threshold can also be computed with a multiplication.
Moreover, I believe that the right implementation should use inline assembly for the multiply instruction, which provides guaranteed performance. The "nearlydivisionless" function avoids the inline assembly by writing some excessively complicated expression, with the hope that a clever compiler will optimize it by simplifying it to something equivalent with the correct assembly instruction. I do not believe that relying on very clever compilers is the right choice for achieving the portability of a program. On the contrary, even when recompiling the program for the same ISA, a new version of the same compiler may not do the desired optimization, causing unpredictable performance problems. If the compiler fails to optimize the 128-bit multiplication written in the "nearlydivisionless" function, by noticing that both 128-bit operands have been recently promoted from 64-bit values, it would do 4 multiplications instead of 1 multiplication.
Like I have said, I think that Daniel Lemire has another more recent article, including better multiplication methods, but I am too lazy to search for it now.
The only error then is when there are no choices or they set the weights to 0, negative or NaN. I am but a lowly musician though. I can't say I always think floats are simple...
I think actually this is an argument for language designers to include a “std.choice” in their standard library that consumes random bytes and correctly performs common ergonomic operations like “get one element at random from this collection”.
(If your standard library tries to make a distinction between “regular random number generators” and “cryptographically secured random number generators”, I think this distinction between “generate random bits” and “make probabilistic choices” is about equally important.)
https://heliosphan.org/zigguratalgorithm/zigguratalgorithm.h...
But, I make games.
https://gegell.github.io/posts/factorio-rng/
My comment was mainly a joke. Since my work is primarily with performance and graphics, I rarely use any proper RNG and usually only need cheap things which look random enough at first glance.
Which NIST SP-800-22 implementation instead of the now-archived paranoid_crypto randomness tests?
paranoid_crypto/docs/randomness_tests.md : https://github.com/google/paranoid_crypto/blob/main/docs/ran...
/? NIST SP-800-22 Rust: https://www.google.com/search?q=NIST+SP-800-22+rust&oq=NIST+...
Sometimes it's possible to whiten random to make it uniform random or normal random;
Whitening transformation: https://en.wikipedia.org/wiki/Whitening_transformation