最近在做伪随机数，关于伪随机数生成的学习笔记整理了一下。其中有很多biased approach可能是之前不会意识到的，学起来还是很带感的。
原文如下，[[]]是双向链接，可以到 我的博客 来获取更好的链接体验 (shameless self plug 😃)
WARNING The author of this routine has been writing random-number generators for many years and has never been known to write one that worked.
Unix Third Edition (V3) manual page
What?
Select a Uniformly Distributed random number from interval [0, s), and this interval may change dynamically.
We will assume that RNG will output a uniformly distributed random number in [0, 2^{32} - 1]
To do it correctly in C++.
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
std::uniform_int_distribution<uint32_t> dist(0, range-1);
return dist(rng);
}
Biased Approach
Classic Modulo
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
return rng() % range;
}
- Biased because
rand() produces numbers in the range [0..2^32), for numbers that can not perfectly divide 2^32, some of them have a lower change to be selected.
- Eg: 52 does not perfectly divide
2^32, and 2^32 mod 52 = 48, so 49..52 have a lower change to be selected.
- One other concern we might have is that some generators have weak low-order bits. That is, the lower bit does not pass the statistical tests (for every bit it must appear 50/50). When we perform
% 52, we are essentially passing the lowest bit straight through the output.
Floating Point Multiply
static uint32_t bounded_rand(rng_t& rng, uint32_t range) {
double zeroone = 0x1.0p-32 * rng();
return range * zeroone;
}
- Note that
0x1.0p-32 is a binary floating point constant for 2^{-32}
- This approach is just as biased as the classic modulo approach, but the bias manifests itself differently. For example, if we were choosing numbers in the range
[0..52), the numbers 0, 13, 26 and 39 would appear once less often than the others.
- That is, it does not uniformly map input to an even output.
Int Multiply
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
uint32_t x = rng();
uint64_t m = uint64_t(x) * uint64_t(range);
return m >> 32;
}
- It might appear that this version requires 64-bit arithmetic, but on x86 CPUs, a good compiler will compile this code to a 32-bit
mult instruction (which yields two 32-bit outputs, one of which is the return value). We should expect this version to be fast, but biased in exactly the same way as the floating-point multiplication method.
- Biased for the exact same reason as the floating-point multiplication method.
Unbiased
The most common method to achieve fairness is through rejection.
Division with Rejection
Instead of calculating x * range / 2**32, we calculate x / (2**32 / range), then reject the uneven condition, essentially reject the last skewed biased (last 48 values).
For example, in the case of drawing a card using a 32-bit PRNG, we can generate a 32-bit number and divide it by 2^32/52 = 82,595,524 to choose our card. This method works when the random value from the 32-bit PRNG is less than 52 × 82595524 = 2^32/32 – 48. Since division do not have remainder.
If the random value from the PRNG is one of the final 48 values at the top end of the range of the generator, it must be rejected and another sought.
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
// calculates divisor = 2**32 / range
uint32_t divisor = ((-range) / range) + 1;
if (divisor == 0) // overflow, it's really 2**32
return 0;
for (;;) {
uint32_t val = rng() / divisor;
if (val < range)
return val;
}
}
- Note that range is unsigned. For an unsigned integer, the negation of a number n is 2^{32} - n.
- Dividing this by the range gives a value just less than 2^{32} / \text{range}, so we add 1 to it to get the correct divisor.
Debiased Modulo (Twice)
- Reject the first 48 value (that are biased)
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
// calculates 2**32 % range
uint32_t t = (-range) % range;
for (;;) {
uint32_t r = rng();
if (r >= t)
return r % range;
}
}
The code for OpenBSD's arc4random_uniform (which is also used on OS X and iOS) adopts this strategy.
Debiased Modulo (Once)
Java adopts a different approach to generating a number in a range that manages to use only one modulo operation except in the relatively rare occasions when rejection occurs. The code is
static uint32_t bounded_rand(rng_t& rng, uint32_t range) {
uint32_t x, r;
do {
x = rng();
r = x % range;
} while (x - r > (-range));
return r;
}
It may take a little thought to figure out why this variant works. Unlike the previous version based on remainders, which removes bias by removing some of the lowest values from the underlying generation engine, this version first uses x - r to make it multiple of range, then filter out value that lands in 2**32 - range as before to make it unbiased.
Debiased Integer Multiplication — Lemire's Method
Much as we removed bias from the modulo technique, we can also remove bias from the integer multiplication method. This method was devised by Lemire.
The highlight is that it does not require any modulo arithmetic, which makes it use less CPU cycles and should be fast enough. This is also the algorithm used by Go's rand package.
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
uint32_t t = (-range) % range;
do {
uint32_t x = rng();
uint64_t m = uint64_t(x) * uint64_t(range);
uint32_t l = uint32_t(m);
} while (l < t);
return m >> 32;
}
Bitmask with Rejection (Unbiased) — Apple's Method
Our final approach avoids division and remainder operations entirely. Instead, it uses a simple masking operation to get a random number in the range [0..2^k) where k is the smallest value such that 2^k is greater than the range. If the value is too large for our range, we discard and try another. Code below:
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
uint32_t mask = ~uint32_t(0); // all bits to 1
--range; // Decrement the range by 1 because the range [0..range-1] needs to be inclusive.
// clz: count leading zero
// range|1 so we have at least 1 bit mask
// right shift effectively creating a mask that matches the bit length of range.
mask >>= __builtin_clz(range|1);
uint32_t x;
do {
// since rng() return bits uniformly, x should also have uniformly distributed bits
x = rng() & mask;
} while (x > range);
return x;
}
This approach was adopted by Apple when (in the macOS Sierra release) they made their own revision to the code for arc4random_uniform.
Use Cases
![[Knuth Shuffle]]
- [[Random Sampling]]
- [[Mersenne Twister]]
PRNG VS RNG
- Pseudo: the input space (seed) is limited. Thus, we can predict the result.
- CPRNG: crypto-strength version of
PRNG. It takes a 256-bit seed rather than a 32-bit seed.
Links