Fast way to generate pseudo-random bits with a given probability of 0 or 1 for each bit

user_5404527 0
Viewed 6316

Normally, a random number generator returns a stream of bits for which the probability to observe a 0 or a 1 in each position is equal (i.e. 50%). Let's call this an unbiased PRNG.

I need to generate a string of pseudo-random bits with the following property: the probability to see a 1 in each position is p (i.e. the probability to see a 0 is 1-p). The parameter p is a real number between 0 and 1; in my problem it happens that it has a resolution of 0.5%, i.e. it can take the values 0%, 0.5%, 1%, 1.5%, ..., 99.5%, 100%.

Note that p is a probability and not an exact fraction. The actual number of bits set to 1 in a stream of n bits must follow the binomial distribution B(n, p).

There is a naive method that can use an unbiased PRNG to generate the value of each bit (pseudocode):

generate_biased_stream(n, p):
  result = []
  for i in 1 to n:
    if random_uniform(0, 1) < p:
      result.append(1)
    else:
      result.append(0)
  return result

Such an implementation is much slower than one generating an unbiased stream, since it calls the random number generator function once per each bit; while an unbiased stream generator calls it once per word size (e.g. it can generate 32 or 64 random bits with a single call).

I want a faster implementation, even it it sacrifices randomness slightly. An idea that comes to mind is to precompute a lookup table: for each of the 200 possible values of p, compute C 8-bit values using the slower algorithm and save them in a table. Then the fast algorithm would just pick one of these at random to generate 8 skewed bits.

A back of the envelope calculation to see how much memory is needed: C should be at least 256 (the number of possible 8-bit values), probably more to avoid sampling effects; let's say 1024. Maybe the number should vary depending on p, but let's keep it simple and say the average is 1024. Since there are 200 values of p => total memory usage is 200 KB. This is not bad, and might fit in the L2 cache (256 KB). I still need to evaluate it to see if there are sampling effects that introduce biases, in which case C will have to be increased.

A deficiency of this solution is that it can generate only 8 bits at once, even that with a lot of work, while an unbiased PRNG can generate 64 at once with just a few arithmetic instructions.

I would like to know if there is a faster method, based on bit operations instead of lookup tables. For example modifying the random number generation code directly to introduce a bias for each bit. This would achieve the same performance as an unbiased PRNG.

Edit March 5

Thank you all for your suggestions, I got a lot of interesting ideas and suggestions. Here are the top ones:

  • Change the problem requirements so that p has a resolution of 1/256 instead of 1/200. This allows using bits more efficiently, and also gives more opportunities for optimization. I think I can make this change.
  • Use arithmetic coding to efficiently consume bits from an unbiased generator. With the above change of resolution this becomes much easier.
  • A few people suggested that PRNGs are very fast, thus using arithmetic coding might actually make the code slower due to the introduced overhead. Instead I should always consume the worst-case number of bits and optimize that code. See the benchmarks below.
  • @rici suggested using SIMD. This is a nice idea, which works only if we always consume a fixed number of bits.

Benchmarks (without arithmetic decoding)

Note: as many of you have suggested, I changed the resolution from 1/200 to 1/256.

I wrote several implementations of the naive method that simply takes 8 random unbiased bits and generates 1 biased bit:

  • Without SIMD
  • With SIMD using the Agner Fog's vectorclass library, as suggested by @rici
  • With SIMD using intrinsics

I use two unbiased pseudo-random number generators:

I also measure the speed of the unbiased PRNG for comparison. Here are the results:


RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 16.081 16.125 16.093 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 0.778 0.783 0.812 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.176 2.184 2.145 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.129 2.151 2.183 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

SIMD increases performance by a factor of 3 compared to the scalar method. It is 8 times slower than the unbiased generator, as expected.

The fastest biased generator achieves 2.1 Gb/s.


RNG: xorshift128plus

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.300 21.486 21.483 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 22.660 22.661 24.662 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 1.065 1.102 1.078 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 4.972 4.971 4.970 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 4.955 4.971 4.971 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical   : 104,857,600

For xorshift, SIMD increases performance by a factor of 5 compared to the scalar method. It is 4 times slower than the unbiased generator. Note that this is a scalar implementation of xorshift.

The fastest biased generator achieves 4.9 Gb/s.


RNG: xorshift128plus_avx2

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.754 21.494 21.878 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 54.126 54.071 54.145 [Gb/s]
Number of ones: 536,874,540 536,880,718 536,891,316
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 1.093 1.103 1.063 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 19.567 19.578 19.555 [Gb/s]
Number of ones: 104,836,115 104,846,215 104,835,129
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 19.551 19.589 19.557 [Gb/s]
Number of ones: 104,831,396 104,837,429 104,851,100
Theoretical   : 104,857,600

This implementation uses AVX2 to run 4 unbiased xorshift generators in parallel.

The fastest biased generator achieves 19.5 Gb/s.

Benchmarks for arithmetic decoding

Simple tests show that the arithmetic decoding code is the bottleneck, not the PRNG. So I am only benchmarking the most expensive PRNG.


RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)

Method: Arithmetic decoding (floating point)
Gbps/s: 0.068 0.068 0.069 [Gb/s]
Number of ones: 10,235,580 10,235,580 10,235,580
Theoretical   : 10,240,000

Method: Arithmetic decoding (fixed point)
Gbps/s: 0.263 0.263 0.263 [Gb/s]
Number of ones: 10,239,367 10,239,367 10,239,367
Theoretical   : 10,240,000

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 12.687 12.686 12.684 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 14.536 14.536 14.536 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 0.754 0.754 0.754 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.094 2.095 2.094 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.094 2.094 2.095 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

The simple fixed point method achieves 0.25 Gb/s, while the naive scalar method is 3x faster, and the naive SIMD method is 8x faster. There might be ways to optimize and/or parallelize the arithmetic decoding method further, but due to its complexity I have decided to stop here and choose the naive SIMD implementation.

Thank you all for the help.

bit-manipulation random c algorithm c++
11 Answers

Although this question is 5 years old, I believe I have something of value to add. While SIMD and arithmetic decoding are undoubtedly great techniques, its hard to ignore that the bitwise boolean arithmetic suggested by @mindriot is very simple and easy to grasp.

However, it's not immediately apparent how you would go about efficiently and quickly implementing this solution. For 256 bits (0.00390625) of resolution, you could write a switch statement with 256 cases and then manually determine the required boolean expression by hand for each case. It would take a while to program this but it will compile down to a very fast jump table in C/C++.

But, what if you want 2^16 bits of resolution, or even 2^64? The latter is a resolution of 5.4210109E-20, more precise than most of us would ever need. The task is absolutely impossible by hand, but we can actually construct a small virtual machine to do this quickly in just 30 lines of C code.

Let's construct the machine for 256 bits of resolution. I'll define probability = resolution/256. e.g., when resolution = 64, then probability = 0.25. As it turns out, the numerator (resolution) actually implicitly encodes the required boolean operations in its binary representation.

For example, what expression generates probability = 0.69140625 = 177/256? The resolution is 177, which in binary is 10110001. Let AND = 0 and OR = 1. We start after the first nonzero least significant bit and read toward the most significant bit. Map the 0/1 to AND/OR. Thus, starting from b1 and reading right to left, we generate the boolean expression (((((((b1 and b2) and b3) and b4) or b5) or b6) and b7) or b8). A computer-generated truth table will confirm 177 cases yield True. To give another example, probability = 0.4375 = 112/256 gives the resolution in binary as 01110000. Reading the 3 bits in order after the first non zero LSB (011) gives ((b1 | b2) | b3) & b4).

Since all we need are the two AND and OR operations, and since the resolution encodes the exact boolean expression we need, a virtual machine can be programmed which interprets the resolution as bitcode. AND and OR are just opcodes that act immediately on the output of an unbiased random number generator. Here is my sample C code:

uint64_t rng_bias (uint64_t *state, const uint8_t resolution)
{
    if (state == NULL) return 0;
    
    //registers
    uint64_t R0 = 0;
    uint8_t  PC = __builtin_ctz(resolution|0x80);
    
    //opcodes
    enum
    {
        OP_ANDI = 0,
        OP_ORI  = 1,
    };
    
    //execute instructions in sequence from LSB -> MSB
    while (PC != (uint8_t) 0x8)
    {
        switch((resolution >> PC++) & (uint8_t) 0x1)
        {
            case OP_ANDI:
                R0 &= rng_generator(state);
                break;
                
            case OP_ORI:
                R0 |= rng_generator(state);
                break;
        }
    }
    
    return R0;
}

The virtual machine is nothing more than 2 registers and 2 opcodes. I am using GCC's builtin function ctz which counts the trailing zero bits so that I can easily find the first nonzero LSB. I bitwise-or the ctz argument with 0x80 because passing zero is undefined. Any other decent compiler should have a similar function. Notice, that unlike the examples I showed by hand, the VM interprets the bitcode starting on the first nonzero LSB, not after. This is because I need to make at least one call to the PRNG to generate the base p=0.5 and p=0.0 cases.

The state pointer and rng_generator() calls are used to interface with your random number generator. For example, for demonstration purposes I can use Marsaglia's Xorshift64:

uint64_t rng_generator(uint64_t *state)
{
    uint64_t x = *state;
    
    x ^= x << 13;
    x ^= x >> 7;
    x ^= x << 17;
    
    return *state = x;
}

All the user/you need to do is manage a separate uint64_t state variable, which must be appropriately seeded prior to using either function.

It is extremely easy to scale to 2^64 bits or whatever other arbitrary resolution desired. use ctzll instead for unsigned long long arguments, change the uint8_t types to uint64_t, and change the while loop check to 64 instead of 8. That's it! Now with at most 64 calls to the PRNG, which is fairly fast, we have access to 5.4210109E-20 resolution.

The key here is that we get the bitcode practically for free. No lexing, parsing, or any other typical VM interpreter tasks. The user provides it via the resolution, without ever realizing. As far as they're concerned, its just the numerator of the probability. As far as we, the implementers are concerned, its nothing more than a string of bitcode for our VM to interpret.

Explaining why the bitcode works requires a whole different and much longer essay. In probability theory, the problem is to determine the generating event (the set of all sample points) of a given probability. Not unlike the usual inverse CDF problem for generating random numbers from a density function. In a computer science viewpoint, in the 256 bit resolution case, we are traversing a depth-8 binary tree where each node represents a probability. The parent node is p=0.5. Left traversal indicates AND operations, right traversal indicates OR. The traversal and node depth maps directly to the LSB->MSB bit encoding that we discussed several paragraphs before.

user_84226480
Related